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ABSTRACT 


This report describes the development of new numerical methods and new constitutive 
models for the simulation of hypervelocity impact effects on spacecraft. The research has 
included parallel implementation of the numerical methods and material models developed 
under the project. Validation work has included both one dimensional simulations, for com- 
parison with exact solutions, and three dimensional simulations of published hypervelocity 
impact experiments. The validated formulations have been applied to simulate impact ef- 
fects in a velocity and kinetic energy regime outside the capabilities of current experimental 
methods. The research results presented here allow for the expanded use of numerical simu- 
lation, as a complement to experimental work, in future design of spacecraft for hypervelocity 
impact effects. 
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INTRODUCTION 


The seven chapters which follow describe numerical methods and constitutive modeling 
work which offers improved capabilities for the simulation of hypervelocity impact effects. 

• Chapter 1 describes the development of a new hybrid-particle finite element method 
for hypervelocity impact simulation. 

• Chapter 2 describes validation of the hybrid-particle finite element method developed 
in the first chapter, in simulations of published hypervelocity impact experiments. 

• Chapter 3 describes validation of the hybrid-particle finite element method developed 
in the first chapter, in simulations of rigid body dynamics problems with known exact 
solutions. 

• Chapter 4 describes an improved formulation of the hybrid-particle finite element 
method developed in the first chapter, simplifying the method and improving its com- 
putational efficiency. 

• Chapter 5 describes the development of a new rate dependent elastic-plastic model for 
Kevlar, and its application in the simulation of hypervelocity impact effects on com- 
posite orbital debris shielding, like that deployed on the International Space Station. 

• Chapter 6 describes the development of a new rate dependent anisotropic elastic-plastic 
model for reinforced carbon-carbon, and its application in the simulation of hyperve- 
locity impact effects on the wing leading edge of the Space Shuttle. 

• Chapter 7 describes the application of hybrid particle-finite element methods in the 
simulation of foam impact effects on the Space Shuttle thermal protection system. 

The Conclusions section suggests opportunities for future research, and Appendix 1 tab- 
ulates some specific simulation data for reinforced carbon-carbon. 
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AN ELLIPSOIDAL PARTICLE-FINITE ELEMENT METHOD 


FOR HYPERVELOCITY IMPACT SIMULATION 

Ravishankar Shivarama 3 and Eric P. Fahrenthold 4 

Department of Mechanical Engineering, 1 University Station C2200 
University of Texas, Austin, TX 78712, USA 

A number of coupled particle- element and hybrid particle- element methods have been de- 
veloped for the simulation of hypervelocity impact problems, to avoid certain disadvantages 
associated with the use of pure continuum based or pure particle based methods. To date 
these methods have employed spherical particles. In recent work a hybrid formulation has 
been extended to the ellipsoidal particle case. A model formulation approach based on La- 
grange’s equations, with particles entropies serving as generalized coordinates, avoids the 
angular momentum conservation problems which have been reported with ellipsoidal smooth 
particle hydrodynamics models. 

KEYWORDS: particle methods, finite element methods, impact simulation 

INTRODUCTION 

A review of the literature on hypervelocity impact simulation shows that most research in 
this field has focused on the development and application of continuum hydrocodes, using ei- 
ther Eulerian, Lagrangian, or Arbitrary Lagrangian-Eulerian (ALE) formulations [1,2,3]. In 
the last decade attention has shifted towards particle based methods [4], in particular smooth 
particle hydrodynamics (SPH). Although both continuum and particle based methods have 
demonstrated excellent results in a range of practical applications, these methods are not 
without problems. As a result some recent research has formulated coupled particle-element 
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[5] or hybrid particle-element [6] methods, aimed at avoiding certain disadvantages of the use 
of pure continuum based methods or pure particle based methods in hypervelocity impact ap- 
plications [7]. The present paper develops a generalized hybrid method, extending the work 
of Fahrenthold and Horban [6] , presenting for the first time a combined particle-element for- 
mulation based on a nonspherical particle geometry. Nonspherical particle formulations are 
also of interest in certain molecular dynamics [8] and astrophysics [9,10] applications, where 
material or flowheld geometry has suggested the use of repulsion potentials or interpolation 
kernels with ellipsoidal forms. 

Most previous hypervelocity impact simulation work has employed Eulerian finite volume 
and Lagrangian finite element methods [1,11] and more recent extensions of these methods 
to Arbitrary Lagrangian- Eulerian frames [12]. Although continuum hydrocodes are accurate 
and efficient simulation tools for many important problems, they are not well suited for use in 
all hypervelocity impact applications. Eulerian codes can accommodate arbitrarily general 
contact-impact, but their approximate models of material strength effects are best suited to 
a very high velocity impact regime. Lagrangian codes incorporate very accurate models of 
material strength effects, but their slideline based contact-impact algorithms are best suited 
to a relatively low velocity impact regime. ALE methods allow for intelligent tailoring of the 
mesh, ideally avoiding (but not eliminating) the preceding limitations. 

None of the aforementioned continuum formulations is ideally suited to applications where 
the generation, transport, and contact-impact dynamics of material fragments is of central 
concern. When material fragments become much smaller than the finite volume cell size, 
Eulerian flow held interpolations may not accurately represent the physical state of the highly 
comminuted medium. Hence the adaptive introduction of a high resolution Eulerian mesh 
is required. Similarly the use of Lagrangian finite elements to represent a highly fragmented 
medium calls for the introduction of a very extensive and adaptive slideline mesh. It is 
not surprising that these continuum formulations can be difficult to apply in cases where 
fragmentation dynamics are of central interest, given the highly discontinuous nature of the 
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field variables in such systems. 

Recognizing the disadvantages of continuum formulations in certain hypervelocity im- 
pact applications, a number of different particle methods [13,14,15] have been introduced. 
In general these methods replace or augment the conventional continuum kinematics with 
particle based interpolations, to more efficiently represent fragment generation, transport, 
and contact-impact effects in applications where modeling of the latter physics is of central 
concern. The overwhelming majority of particle based hypervelocity impact models have 
employed SPH techniques and spherical kernels. A notable exception is the work of Shapiro 
and co-workers [9,10], who extended the basic SPH method to ellipsoidal kernels. However 
those authors indicate that their method fails to conserve angular momentum, and they of- 
fer no solutions to the tensile instability, numerical fracture, and other problems which have 
hindered the effective use of SPH methods in some applications. The formulation presented 
here does not involve element-to-particle conversion, as discussed for example by Johnson et 
al. [14], who converted distorted finite elements into SPH particles. The latter conversion, 
if it occurs before failure of the element, may introduce the aforementioned SPH tensile 
instability and and numerical fracture problems. 

The present paper formulates a hybrid ellipsoidal particle-finite element model for hyper- 
velocity impact simulation. Unlike the coupled particle-element formulations of Hayhurst et 
al. [16] and other authors, which use particles to model some structures and finite elements 
to model other structures, the present paper employs particles and elements everywhere. The 
simultaneous use of particles and elements is not redundant, since they are used to represent 
distinct physics. The particles model all inertia and contact-impact effects, and are associ- 
ated with interpolation functions which represent compressed states. The elements model all 
strength effects, namely tension and elastic-plastic shear. This approach has several advan- 
tages. The tensile instability and numerical fracture problems common to SPH formulations 
are avoided, since interpolation kernels are not used to represent interparticle tension forces, 
and since material strength effects give rise to interparticle forces only among reference con- 
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figuration neighbors. Although true Lagrangian kinematics are used to calculate tensile and 
shear forces, the use of particles in all inertia and contact-impact calculations means that 
material failure is simply accommodated (without element-to-particle remapping or mass 
and energy discard) via the loss of element cohesion. Particles not associated with any in- 
tact elements translate and rotate as individual fragments, in response to contact-impact 
loads. No slidelines are used, so that contact-impact of all intact and fragmented material 
is modeled everywhere. Although the authors do not contend that this numerical approach 
is best for all impact simulation problems, it offers important advantages in hypervelocity 
impact applications. In the latter field an interest in multi-structure perforation, erosion, 
fragmentation, and very general contact-impact modeling is not unusual. 

The model formulation work described in the sections which follow is different from that 
normally employed to develop hydrocode models of either the continuum or particle type. 
Two particular features should be mentioned. First an energy method (Lagrange’s equations) 
is used to develop the semi-discrete model, so that no reference is made to any partial 
differential equations. The introduction of entropy variables as generalized coordinates makes 
it possible to account for very general thermal dynamics while using a discrete Lagrangian 
approach. Secondly the formulation introduces Euler parameter states [17] to account for 
the rotational motion of the nonspherical particles. An Euler parameter description of the 
particle kinematics is desired in order to avoid the singularities of Euler angle based models. 
Although the use of Euler parameters leads nominally to a differential-algebraic formulation, 
the development which follows introduces momentum state variables, to reduce the order 
of the system and thereby eliminate the Lagrange multipliers associated with the Euler 
parameter constraint. The model formulation work described here therefore generalizes in 
a geometric sense the particle-element work of Fahrenthold and Horban [6] and generalizes 
in a thermomechanical sense the Euler parameter based mechanical models of Chang and 
Chou [18] and others [19,20]. 

The organization of this paper may be summarized as follows. The second and third 
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sections define the particle and element kinematics and the interpolations used for all held 
variables. The fourth and fifth sections develop kinetic energy and thermomechanical po- 
tential energy functions for the system. The sixth section describes the dissipative process 
models, for plastic formation and damage evolution. The seventh and eighth sections present 
artificial viscosity and artificial heat diffusion models, similar to those used in other shock 
physics codes, and develop entropy evolution equations which take the place of the energy 
conservation relations normally employed. The ninth section introduces a virtual work ex- 
pression, to account for externally applied loads. The tenth section combines the canonical 
thermomechanical Lagrange equations with the multiple nonholonomic constraints devel- 
oped in the preceding sections, and derives an unconstrained explicit state space model 
for the particle-element system. The final section presents two example problems, show- 
ing good agreement of simulations performed using the model developed here to published 
data for three dimensional impact experiments. Additional validation work, for simulations 
performed using a spherical kernel, are described by Fahrenthold and Shivarama [21], who 
also provide details on the numerical implementation and test results showing numerical 
convergence and good parallel speedup. 


PARTICLE KINEMATICS 


The inertia distribution in the modeled system is represented by a collection of n ellipsoidal 
particles, with the mass of the Ah particle, whose position and orientation are described 
by a center of mass position vector (c^O) and an Euler parameter vector (e^) with component 
forms 


(0 - r Jd AO AO i t 


c = c 


(0 = [ e (0 ao AO AO } T 


e K ' = 


o c i 


( 1 ) 


where T denotes the transpose. The components of the center of mass position vector 
are described in a fixed Cartesian coordinate system. The Euler parameters [17] define an 
orthogonal rotation matrix (R(0) for each particle 


R« = hW g wt 


( 2 ) 
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which relates the components of any vector (a), described in a body fixed co-rotating coor- 
dinate system, to its corresponding components in the fixed Cartesian system, using 

a = R® a (5) 


The four Euler parameters may be used to compute three Euler angles for each particle, 
and are preferred for use in large rotation dynamics problems, since the Euler angles are 
singular functions. Although the Euler parameters are nonsingular, they must satisfy the 
constraint 

e® T e® = 1 ( 6 ) 

which requires that 

G®e® = 0, G®e® = -G®e®, G® G® T = ! (7) 


where I denotes an identity matrix. 

The Euler parameters and their time derivatives are related to the particle angular velocity 
vectors ( u>® ), with components expressed in the particle’s co-rotating frame, by the well 
known [17] relations 

= 2 G® e®, e® = ^ G® T u>® (8) 

Similarly the anti-symmetric matrix ft, with axial vector u>, which satisfies 

G®v = w (0 x v ( 9 ) 


for any vector v, is related to the Euler parameters and their time derivatives by 

ft ® = 2 G® G® r = -2 G®G® T (10) 
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The next section presents a density interpolation associated with the particles just de- 
scribed, and defines Lagrangian finite elements whose nodal coordinates are components of 
the particle center of mass position vectors. 


DENSITY INTERPOLATION AND FINITE ELEMENTS 


In the material reference configuration, the ellipsoidal particles are arranged in a packing 
scheme defined by a mapping from a body centered cubic unit cell, the mapping composed 
of stretchings in three orthogonal directions, those directions aligned with the principal axes 
of the ellipsoids. The density interpolation for compressed states is 


p ( i) = p f + ^ iy«> (ii) 

j = 1 


where is the continuum density at the centroid of the fth particle, pd ! is the reference 
density for the ith particle, and is the interpolation kernel 




1 
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( 12 ) 




(13) 


The symbols S tJ and u denote respectively the Dirac delta function and the unit step function 


u(x ) 


0, x < 0 

1, x > 0 


(14) 


The constants a and (3 are effective separation distances, in a unit cell, measured respectively 
between body centered particles and between a body centered particle and a particle located 
at a cell vertex 


a = 




(15) 


Note that the lead coefficient in the expression for is due to the presence of eight 

nearest neighbors in the reference configuration. 


The function is an ellipsoidal coordinate, defined in a frame which co-rotates with 


the jth particle 
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where h\ 3 \ h^\ h\p are the half-lengths of the principal axes of the j'th particle. 

Several properties of the density interpolation and kernel functions just described should 
be noted. The interpolation provides an exact Lagrangian description of the variation of 
density under uniform compression, and incorporates a lower bound equal to the particle ref- 
erence density, so that tensile instabilities are avoided. The kernel is singular, so that no two 
particle centers overlap, avoiding particle streaming effects. In the reference configuration 
neighbor particles do not contribute to the density calculation, so that no special treatment 
of surface particles with incomplete neighbor sets is required. The density dependence of 
the compact kernel support means that particles interact only with near neighbors, and that 
mechanical interaction with remote neighbors through intervening matter is avoided. Finally 
note that the dimensionless kernel used here is closer in functional form and physical inter- 
pretation to the repulsion potentials used in molecular dynamics, than to the dimensional 
kernels with cubic spline form used in most SPH formulations. 

The center of mass coordinates of vertex centered particles are nodal coordinates for hex- 
ahedral finite elements, which are used to account for tension and elastic-plastic shear in 
cohesive solid materials. These elements are used with one point integration, as described 
by Hallquist [22], although in the present case the element packing scheme and density in- 
terpolation preclude the development of hourglass modes. Since large deformation problems 
are of interest in hypervelocity impact applications, the element nodal coordinates are used 
here to compute a Jacobian (J^) and a Lagrangian deviatoric strain tensor (E^) for each 


9 


element, defined by [23] 


Jti) = det F (j) , E U) = ^ ( C (i) - 1 ) (18) 

where F (y ^ is the deformation gradient for the jth element and 

C°) = F°' )T F 0) , F {j) = ( det F F W) (19) 

The plasticity model discussed in a later section assumes an additive decomposition of the 
deviatoric strain, with the elastic deviatoric strain tensor (E ft W) f or the jth element dehned 

by 

E e(j) = E 0) - E P(J ') (20) 

so that the evolution equations for the plastic strain tensor (E p bl) must satisfy the isochoric 
plastic deformation constraint 

tr (c p( b- T C p(j) ^ = 0, E p(i ) = I (cpCj): _ i) (21) 

It should be emphasized that the numerical modeling methodology developed in this paper 
applies for a wide range of element types and plasticity models. The particular element 
kinematics and constitutive equations used in this paper are representative formulations 
which account for large deformation kinematics. 

KINETIC ENERGY 

The kinetic co-energy for the system is 

n 

T* = ^2 T * (i) ( 22 ) 

1 = 1 

where T*W j s the kinetic co-energy for the ith particle 

T*(*) = - m {i) c w T c w + - T jW w (i) (23) 

2 2 v ’ 

and the components of the inertia matrix Jb) are calculated in a body fixed frame. Since 

the components of the angular velocity vector are quasi-velocities, equations (8) are used to 
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rewrite the second term 




- c ^ T c« 
2 


+ 2e (i)T G (i)T J (i) G w e w 


(24) 


which identifies the center of mass coordinates and the Enler parameters as generalized 
coordinates in Lagrange’s equations. Note that the relations (7) allow the rotational kinetic 
co-energy to be expressed as 


rp*(i ) 
rot 


2e (i)T G (i) T J (i) G w e w 


(25) 


so that the generalized force associated with the Euler parameter dependence of the kinetic 
co-energy is 

r)T*(h 

k h) = = 4 G« t J«g« e (i) (26) 

<9eW v ’ 

The next section develops a potential energy expression for the particle-element system. 


POTENTIAL ENERGY 

The potential energy of the particle-element system is the sum of the particle internal 
energies and the element stored energies due to tension and elastic shear, and has the general 
form 


^ m(i) u{i) + J 2 ( £/ten(j) jjdev(j) _|_ jjpen(j) j 

3 = 1 


(27) 


i= 1 


where j s an internal energy per unit mass, n e is the number of elements, U ten ^ is a 
stored energy in tension, JJ dev ^' ) is a stored energy in shear, and U pen ^ is a penalty energy 
used to position the body centered particle in each element. The particle internal energy 
density is determined by the mass density and entropy density (s^) 


U® =U {i) (p W ,S W ) , 


s « = 


m 


(0 


(28) 


where S ^ is the total entropy of the fth particle. The functional form of the internal energy 
per unit mass depends on the equation of state, and it is not unusual in such calculations to 
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rely upon tabular data [24]. In any case, the extensive state variables which determine the 
particle internal energy are 

u® =u® (c®,e®,S®) (29) 

The form of the element stored energy functions depends upon the constitutive assump- 
tions. Although the general method developed here admits a wide range of constitutive 
models, simple forms are assumed here, and are described as follows. The stored energy in 
tension is 

U ten U) = I (l - D u) ) k& V} j) (J u) - l) 2 u{ J u) - 1) (30) 

where is a bulk modulus, is the normal damage, and vd 3> is the element reference 
volume. The stored energy in shear is 

U dev{j) = (1 - d ij) ) fi ij) Vj j) tr (E e ^ T E e(j ' } ) (31) 

where ii ( j> is a shear modulus and d {j) is the deviatoric damage. The penalty function is 

U pen{j) = ^ (1 - D( j) ) K® (c (j) - c (i) ) 2 (32) 

where is the center of mass position vector for the body centered particle of the jth 
element, is the average of the center of mass position vectors for the vertex centered 
particles (the element nodes), and the penalty stiffness is 

K {j) = E U) V 0 (j) ^ (33) 

where is Young’s modulus for the element. The damage variables are used to model the 
loss of cohesion due to element failure, and are discussed in the next section. 

The preceding constitutive assumptions combine with the element kinematic relations 

jU) = jU) ( C M) ? E e(j) = E e(j) epOO) (34) 

to yield the following state variable dependence for the system potential energy 

V = V (c w , e w , S (i) , E p(j \ D®, d U) ) (35) 
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The listed arguments, which include the particle entropies, are by definition Lagrangian 
generalized coordinates. 

The system potential energy function defines the conservative generalized forces 


dV 


= g« 


dV 


- = M {i) 

A ±VA > 


dV 


= e (i) 


(36) 


dc® ’ <9eb) ’ dS® 

where 6® is a particle temperature, and gW and are forces and torques which act on 
the ellipsoidal particles. The energy conjugates for the dissipative state variables define a 
deviatoric stress 


as well as the energy release rates 

Y d U) = 


S u) = 


dV 


dV 


y(j) dWW 


pAl) — 


dV 


dD® ’ ddd) 

The next section discusses evolution equations for the dissipative state variables. 


(37) 


(38) 


PLASTICITY AND DAMAGE MODELS 

The evolution equations for the plastic strain tensor and continuum damage variables will 
serve as nonholonomic constraints on the system level Lagrange equations. As discussed in 
the last section, the general formulation developed here admits a wide range of constitutive 
models. The relatively simple dissipative constitutive models assumed here are described 
in this section. The plasticity model is adapted from Fahrenthold and Horban [25], and 
represents the simplest possible accommodation of the isochoric deformation constraint (21). 
The damage evolution relations are adapted from the Eulerian hydrocode work of Silling [26]. 

The flow rule for the plastic strain is 


j. 0) 

E K?) = - - (C p ^) W ^ + W ^ C p ^) 

net) 

W et) — QPU) gptf) _|_ g p(j) QP(J) _ I^ r (QPti) g pU) g p(j) QP{3)^j J 
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where A^) is a scalar multiplier and 

jqO) _ n _|_ r jU) e Pti)\ N< ' j) 


tr (W U)T W U) ) 


1 1/2 


(39) 

(40) 

(41) 
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e pU) = 


i tr (^pU)T-^p(j)^ 


1/2 


( 42 ) 


with rj (j> a strain hardening modulus, a strain hardening exponent, and e p ^ the accu- 
mulated plastic strain. The effective stress tensor which appears in the flow rule is 


gpd) _ ^ _|_ ^d) € pU)\- nW gd) ; 


s “ = ^^) = (1 - d “ )2 ' , “ E, “ (43) 
where the indicated deviatoric stress is power density conjugate to the plastic strain rate, in 
the entropy equality for the solid medium. 

The yield condition is 


yd) _ T d) _ yd) 


T d) _ 


-tr (S p(i)T S p(i) ) 


1 1/2 


( 44 ) 


where r ^' 1 is the second invariant of the effective stress. The yield stress Y is 

y(i) = yd) (1 _ 7 (J)0«(J)) 


( 45 ) 


where is the reference yield stress, 7^ is a thermal softening modulus, and 0 H ^ is the 
maximum historical homologous temperature. The plastic strain increment at each time 
step is determined using a one step iteration procedure [22] with 

-O') _ yd)) fi( T d) _ yd)) 


AA (i) = 


T v 


(1 - dd)) 2 pd) 

The final evolution equations for the plastic strain have the functional form 


( 46 ) 


E rU) _ j;rU) 


( 47 ) 


and are nonholonomic constraints on the particle-element model. 

Hypervelocity impact problems normally involve perforation and fragmentation effects, so 
that the ability to model such physics is essential. I11 the present case normal and deviatoric 
damage variables are introduced, to model the transition from a cohesive solid, characterized 
by intact finite elements, to a comminuted medium, described by the free flow of particles. 
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Free particles interact with cohesive material and with other fragments under general contact- 
impact loads. 

The damage evolutions applied here are 

D® = u( 1 - D (j) ), d ® u( 1 - d (j) ) (48) 

fiAt K h n At v v ' 

where At is the time step, n is the number of time steps used to model the transition from 

an intact to a failed element, and 

A w = max{ u(e p{3) - e p f U) ), u(6^ ax - 9$), u(J [ c j) - J^ n ), u(a & L - } (49) 

The function just defined initiates damage evolution in an element when when the accumu- 
lated plastic strain e p( ^\ maximum temperature 9 max , minimum element Jacobian J ^] n , or 
maximum eigenvalue of the deviatoric stress (Jrnax reach corresponding critical values for the 
failure strain (e^) , melt temprature (6m), maximum compression (jj^), or spall stress 
(aF' ) ). More sophisticated damage evolution relations may be introduced without change 
to the general modeling framework. The damage evolution equations are nonholonomic 
constraints which apply to the system level Lagrange equations. 


ARTIFICIAL VISCOSITY AND HEAT DIFFUSION 


The shock physics problems of interest here call for the introduction of artificial viscosity 
and artificial heat diffusion effects. The standard formulation used here takes the viscosity 
force pb) on the ith particle to be 


pw = ^2 u(i,j) max (°y ij) ) r M , 

3 = i 


r (*d) 


( C W - c«) ) 

I cb) — cb) I 


where the velocity difference ybJ) j s positive for converging particles 

( f.(i) _ f.U) j . r (0) 
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The damping coefficient z/b J ) j s 
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= ^,(0 c WyW3 +p U) 
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where ci l> is the sound speed for the ith particle and the parameters c Q and c\ are dimen- 
sionless linear and quadratic numerical viscosity coefficients. 

The numerical heat diffusion model assumed here takes the thermal power flow for the 
ith particle to be 

n 




1 

2 


QCon(i) 


Y R (i,j) ( 0 (<) ~ 0 {j) ) 

3 = 1 



(53) 

(54) 


where is the specific heat for the it\i particle and the parameter k Q is a dimensionless 
numerical heat diffusion coefficient. 


ENTROPY EVOLUTION EQUATIONS 

The introduction of entropy state variables, which facilitates the use of an energy based 
modeling approach, means that the conservation of energy relations normally incorporated 
in shock physics models are replaced here by evolution equations for the particle entropies. 
These evolution equations are 

i ) nirr(i) 


gcon(i) 


(55) 


where the irreversible entropy production (S'* rr M) is due to plastic deformation, damage 
evolution, and viscous disspation 




f»( i ) T c(*) 


+ E^‘ 


<j) 


+ r + V„ w tr (s u)T E pW 


(56) 


3 = 1 J 

and is the fraction of the dissipation in element j associated with the ith particle. The 
conduction entropy flow ( 1 S' con (*)) is due to artificial heat diffusion 

n gti) 

gcon(i) = 9 (i)-l QCon(i) = ^ ( 1 - — ) (57) 

3 = 1 

Since the dissipated energy is not lost but rather transduced to the thermal domian, the en- 
tropy evolution equations are nonholonomic constraints in the thermomechanical Lagrangian 
formulation. 
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VIRTUAL WORK 


Although external loads are normally not considered in hypervelocity impact applications, 
for completeness this section develops a virtual work expression. The virtual work for the 
system is a summation over the particles 

n 

sw = 6w(i) (58) 

i = 1 

In this case the particle angular velocities are the time derivatives of rotational quasi- 
coordinates (qW) 

q W = (59) 

so that the virtual work for the ith particle is 

5W {i) = f (i)T hc w + T (t)T hq w (60) 

where and TO are externally applied forces and torques. In view of equation (8), virtual 
changes in the quasi-coordinates and the Euler parameters are related by 

hq (i) = 2 G w he w (61) 

so that the particle virtual work expression is 

5W (i) = f 1 ' i)T 5c {i) + 2 [ G wt T w ] T he w (62) 

The coefficients of the virtual changes in the generalized coordinates are by definition gen- 
eralized nonconservative forces in the system level Lagrange equations. 

LAGRANGE’S EQUATIONS 

The stored energy functions, constraint equations, and virtual work expression developed 
in the preceding sections may be combined with the canonical Lagrange equations, to obtain 
an ODE model for the particle-element system. 


17 


The canonical Lagrange equations are 


d / dT* \ dT* 


dt 

\dc® J 

dc® 

d , 

( dT* \ 

d T* 

dt ' 

yde® J 

de® 

™ -0.(0 
dS® ^ : 

i 

9V — Q d ® 

3d® ^ ’ 


dc® 

= Q c ® 

(63) 

dV 

de® 

= Q e ® 

(64) 

dV 
<9E pCj) 

= QpCj) 

(65) 

dV 

dD® 

= Q D ® 

(66) 


where Q c ^, Q e ^, Q s ®, Q d ®, Q D ^\ and Q p ® are generalized forces determined by the 
constraints and the virtual work. Note that the rate form of the Euler parameter constraint 
for the it\i particle is 

e (i)T e (i) = 0 (67) 

Introducing a Lagrange multiplier j e ® for each Euler parameter constraint as well as La- 
grange multipliers 'y d ® j ^d{i) ^ anc [ yvd) f or the entropy, normal damage, deviatoric 
damage, and plastic strain evolutions equations, the generalized nonconservative forces are 
found to be 


Q c ® 

_ f (i) / rv(i) 

6® 


(68) 

Q e ® 

= 2 G^ r T® + y 6 

(*) e (») 

(69) 

Q s ® 
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If these results and the degenerate Lagrange equations (65) and (66) are combined with the 
previously derived expression for the kinetic co-energy, the momentum balance relations take 
the form 


d 

dt 


( m ® c^) + 


dV 

<9cb) 


f (*) _ p(») 


(74) 


d 

dt 


(4G (i)r J w G w e (i) ) -k w + 


dV 

<9eh) 


2 QbA'jM ^e(i) e (i) 


Introducing the momentum variables 


p ('0 = m (0 c w 


h w = J w 


(75) 


(76) 


and premultiplying the angular momentum balance expression (75) by | G^, which [using 
equation (7)] removes the last term, yields the final state space formulation 


pW 

_ _g(») _ p(») _|_ f(<) 

(77) 

hW 

= -f2 (i) h (i) - -G (i) M (i) + T (i) 
2 

(78) 

c« 

= m (<)_1 p (<) 

(79) 

e« 

— -G( i ) T j(*) -1 h( i ) 
2 

(80) 


Augmented by the evolution equations for S^'\ D (j \ d! J \ and E p b), these are explicit first 
order equations for the particle-element system. With appropriate initial conditions they 
may be integrated to obtain solutions for a wide range of hypervelocity impact problems. 


EXAMPLE SIMULATIONS 

This section compares the results of simulations performed using a numerical implementa- 
tion [27] of the model developed here, to published experimental results for two hypervelocity 
impact problems. The simulations employed a Mie-Grunsisen equation of state and 

c 0 = 0.001, ci = 0.0, k Q = 0.1, e p f = 1.0 (81) 
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with material properties taken from Steinberg [28]. 

The first example involves the impact of a uranium alloy rod on a steel plate, 0.64 cm in 
thickness, in an experiment described by Hertel [29]. The cylindrical projectile is 0.767 cm 
in diameter, with a length to diameter ratio of ten, and impacts the plate at an obliquity of 
73.5 degrees and a velocity of 1.21 kilometers per second. The velocity of the target plate is 
0.217 kilometers per second. The simulation employed 712,929 particles, with the projectile 
composed of spherical particles and the target plate composed of ellipsoidal particles with 
aspect ratios of 1.5:1. 5:1.0. Figure 1 shows an element plot of the initial configuration, while 
Figure 2 shows a particle plot of the simulation results at 100 microseconds after impact. 
The simulation results for the rod erosion (27.5 percent) and residual rod velocity (1.10 
km/s) show good agreement with the corresponding experimental values (27.6 percent and 
1.07 km/s). 

The second example involves the impact of a tungsten rod on a steel plate, 0.95 cm in 
thickness, in an experiment described by Yatteau et al. [30]. The cylindrical projectile is 
0.475 cm in diameter, with a length to diameter ratio of twenty, and impacts the plate at 
an obliquity of 75 degrees and a velocity of 1.83 kilometers per second. The target plate 
is stationary. The simulation employed 671,176 particles, with the projectile composed of 
spherical particles and the target plate composed of ellipsoidal particles with aspect ratios of 
1.5:1. 5:1.0. Figure 3 shows an element plot of the initial configuration, while Figure 4 shows 
a particle plot of the simulation results at 150 microseconds after impact. The simulation 
results for the rod erosion (37.9 percent) and residual rod velocity (1.60 km/s) show good 
agreement with the corresponding experimental values (40 percent and 1.78 km/s). 

In the preceding problems the vast majority of the particles are associated with the target 
plates, so that the use of ellipsoidal particles reduces the total particle count by a factor of 
approximately 2.25, equal to the product of the three indicated aspect ratios. It should be 
noted that the computational advantages of a reduced particle count can in part be offset 
by an increase in the average number of nominal neighbors per particle. This increase is 
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due to the fact that the neighbor search for each particle must consider a spatial volume 
proportional to the cube of the major ellipsoidal axis, despite the fact that a closer approach 
without contact can occur when particle minor axes are appropriately aligned. Experience to 
date indicates that the use of ellipsoidal particles leads to a reduction in memory requirements 
which is proportional to the reduction in particle count, but to a cpu time requirement similar 
to that measured in corresponding simulations with spherical particles. 

CONCLUSION 

The present paper has developed an ellipsoidal particle-finite element method for hyper- 
velocity impact problems, and demonstrated its use in the simulation of three dimensional 
impact experiments. The hybrid particle-element kinematic scheme allows for the use of 
true Lagrangian strength models, while incorporating a completely general description of 
contact-impact effects. Tensile instability, numerical fracture, and angular momentum bal- 
ance problems, often associated with particle formulations, are avoided. Material remaps, 
mass or energy discard, slideline algorithms, and adaptive mesh refinement, often associated 
with continuum models, are also avoided. Derivation of the model is based on a thermome- 
chanical form of Lagrange’s equations, with entropy variables as generalized coordinates. The 
use of an energy based modeling approach facilitates the systematic integration of diverse 
particle based and element based interpolations. Work to date suggests that the method pro- 
vides a numerically robust and accurate approach to the simulation of hypervelocity impact 
problems. 
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Figure 1: Uranium alloy long rod impact on a steel plate, element plot of the initial config- 
uration 
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Figure 2: Uranium alloy long rod impact on a steel plate, particle plot at 100 microseconds 
after impact 
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Figure 3: Tungsten long rod impact on a steel plate, element plot of the initial configuration 
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Figure 4: Tungsten long rod impact on a steel plate, particle plot at 150 microseconds after 
impact 
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Abstract 

The hybrid particle-finite element method of Fahrenthold and Horban [1] , developed for the 
simulation of hypervelocity impact problems, has been extended to include new formulations 
of the particle-clement kinematics, additional constitutive models, and an improved numer- 
ical implementation. The extended formulation has been validated in three dimensional 
simulations of published impact experiments. The test cases demonstrate good agreement 
with experiment, good parallel speedup, and numerical convergence of the simulation re- 
sults. 

Keywords: impact simulation, numerical methods 


1. Introduction 

In previous work Fahrenthold and Horban [1] developed a hybrid particle-finite element 
method for hypervelocity impact simulation, and applied that formulation in the analy- 
sis of several three dimensional problems. The referenced modeling methodolody employs 
Lagrangian finite elements to represent material strength effects, namely tension and elastic- 
plastic shear, and Lagrangian particles to represent inertia, compressed states, and contact- 
impact effects. Particles and elements are used in tandem throughout the simulation, for all 
materials, so that no element-to-particle transformations are required. The introduction of 
both particles and elements is not redundant, since they are employed to account for distinct 
physical effects. A systematic approach to the formulation of this hybrid numerical scheme is 
provided by Hamiltonian mechanics. General thermomechanical dynamics are represented, 
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via the introduction of entropy variables as generalized Hamiltonian coordinates. 

Development of the particle-element method just outlined has been motivated by difficul- 
ties encountered in the application of pure Eulcrian, Lagrangian, or particle based methods 
to the problem of orbital debris shielding design [2]. The latter application calls for simu- 
lation of shock loading and perforation at very high velocities, accurate characterization of 
strength-dependent structural response, efficient modeling of fragment transport, and gen- 
eral descriptions of contact-impact. Recent research suggests that hybrid [3] or coupled [4] 
numerical methods are best suited to simulate orbital debris impact on space structures. 
Numerical methods developed for the latter application may also have application in related 
problems, such as research on the effects of behind armor debris. 


2. Extended Formulation 


The hybrid particle-finite element model of Fahrenthold and Horban [1] has been extended 
to include alternative formulations of the particle-element kinematics, additional constitutive 
models, and an improved numerical implementation. The extended formulation described 
in this section has been evaluated for accuracy, numerical convergence, and parallel speedup 
in a series of three dimensional simulations of published experiments, as described in the 
sections which follow. 

The density interpolation of reference [1] used distinct kernels for reference configura- 
tion nearest neighbors and for all other particles. A simplified alternative is the implicit 
interpolation 
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where p® is the interpolated density, p Q ^ is a reference density, h > ' :i} is a particle radius, r\j 
is a particle center of mass separation distance, the constants a and (3 are determined by the 
body centered cubic particle packing scheme, n is the number of particles, A denotes the unit 
step function, and the summation applies for j ^ i. Note that the preceding interpolation 
yields exact results for hydrostatic compression of a particle set arranged in the selected 
packing scheme, and that in general all neighbor particles interact in the same fashion. 

The particle packing scheme used in the present work leads to eight nearest neighbors 
for body centered particles in the reference configuration. Fahrenthold and Horban [1] used 
the body centered particle associated with each hexahedral element to define subelement do- 
mains. The volumes of these subdomains may be used to determine the interparticle tensile 
forces associated with material dilatation. This six point integration scheme is computa- 
tionally expensive, and leads to hexahedral elements which are stiffer than those employed 
for example in some very successful Lagrangian codes [5] . In the present work the potential 
energy contribution due to tension is written 


n e .. 

U ten = - (1 - D {j) ) {V {j) - l) 2 A (J (j) - 1) + 2 h U) E U) | c U) - c av9U) | 2 } (2) 

3 = 1 
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where n e is the number of elements, V^' 1 is the element reference volume, K/ J) is the element 
bulk modulus, is the element Young’s modulus, is the element Jacobian, D (j> is the 
element normal damage, c-b is the position vector of the body centered particle, and c av9 ^' ) 
is the average position vector of the particles which define the element nodes. This potential 
function represents a one point integration of the element, to quantify tension effects, along 
with a penalty term which positions the body centered particle. Note that interparticle 
tension is (neglecting surface tension) a strength effect, and is therefore associated with 
relative motion of a reference configuration neighbor set (in this case the finite elements) 
and not a time varying neighbor set of the type used to quantify collision forces [6] . 

Shock physics codes can include options for the treatment of energy conservation errors 
which may occur during integration or transformation of the system level model. Examples 
are a default energy discard used in the remap step of CTH [7] and energy errors which 
may arise from contact-impact calculations in ALE [8] or SPH [9,10] formulations. Although 
precise energy conservation can be demonstrated in simple benchmark problems, it can be 
computationally expensive to achieve the same result in the simulation of more complex 
practical problems. For the hybrid particle-element formulation described here, an optional 
energy correction term has been added to the particle entropy evolution equations 


A S err 


e A H err 


3 = 1 
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where A H err is the error in the system Hamiltonian, 6 (j> is a particle temperature, and with 
e = 0.1 the error correction is introduced over ten time steps. The effect of this term is 
to satisfy global energy conservation by introducing an internal energy correction for each 
particle which is proportional to the current particle temperature. 

Simulation results obtained using particle based models are in general influenced by the 
choice of particle packing scheme. The locations and properties of nearest neighbor particles 
define in general a local stiffness distribution for the medium which is anisotropic, as in the 
case of pure crystals [11]. In addition the choice of a particle packing scheme determines 
an effective void ratio for the medium. Hence it is not surprising that particle packing can 
influence simulation results. The preprocessor used in the present work employs a random 
number generator to delete a user specified fraction of body centered particles from the 
original perfect lattice structure. The introduction of these flaws mimics the influence of grain 
orientations, dislocations, and other defects on the mechanical response of real materials. 

The simulations described in references [1] and [3] employed analytic equations of state. 
Extension of the code described in the latter work has included the introduction of interpo- 
lation routines, to accommodate tabulated equations of state in two independent variables. 
Currently the SESAME tables [12] are used for simulations at very high velocities, and are 
accessed in their default form, with pressure and internal energy defined on a density and 
temperature grid. Iteration is therefore used to converge on an internal energy calculated 
from the Gibbs relation 
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where is a particle mass and PW i s a particle pressure. An initial call to the SESAME 
routines is used to establish the reference internal energy associated with the simulation 
initial conditions. 

In addition to the preceding work, significant extensions and applications of the hybrid 
particle-element method discussed here are reported elsewhere. First, Shivarama [13] has 
extended the basic formulation to include ellipsoidal particles, introducing a nonspherical 
kernel, rotational motion of the particles, and Euler parameters as state variables. This 
extension makes possible for example the modeling of thin plate structures at greatly reduced 
particle counts. It should be noted that previous work on nonspherical particle models for 
shock physics problems has been very limited [14,15]. Second, application of the method 
to composite orbital debris shielding problems is reported by Fahrenthold and Park [16], 
including development and numerical implementation of a rate-dependent material model 
for Kevlar. 


3. Validation Simulations 

A series of three dimensional simulations was performed to evaluate the extended formu- 
lation just described. The simulations modeled published experiments conducted at impact 
velocities ranging from one to eleven kilometers per second. Each example problem was 
modeled at two different mesh densities, to investigate numerical convergence of the simu- 
lation results. The selected problems have been previously studied by various investigators, 
to evaluate other numerical methods and computer codes. Tables 1 and 2 provide details on 
the problem parameters and material properties for all the validation simulations. Material 
properties were estimated using data from Steinberg [23]. The first two example problems 
used a Mie-Gruneisen equation of state, while the third example problem used the SESAME 
tables. 


Table 1. Parameters of the example problems 


Parameter 

Sphere 

Long rod 

Debris shield 

Projectile material 

Aluminum 

DU 0.75% Ti 

Aluminum 

Target /shield /wall material 

Aluminum 

Steel 

Aluminum 

First shield thickness (cm) 

na 

na 

0.25 

Second shield thickness (cm) 

na 

na 

0.25 

Target /wall plate thickness (cm) 

0.1143 

0.64 

0.50 

Shield spacing (cm) 

na 

na 

6.0 

Projectile velocity (km/sec) 

6.56 

1.21 

11.0 

Impact obliquity (deg) 

45 

73.5 

45 

Projectile diameter (cm) 

0.953 

0.767 

0.5062 

Projectile length-to-diameter ratio 

na 

10 

4.36 


na = not applicable 
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Fig. 1. Element plot of the initial 
configuration for the sphere impact 
problem. 



Fig. 2. Particle plot of the simulation 
results for the sphere impact problem. 


Table 2. Material properties used in the example problems 


Material property 

Aluminum 

DU 0.75% Ti 

Steel 

Shear modulus (Mbar) 

0.271 

0.74 

0.801 

Reference density (g/cc) 

2.7 

18.62 

7.842 

Initial yield stress (Mbar) 

0.0029 

0.0095 

0.012 

Maximum yield stress (Mbar) 

0.0058 

0.0220 

0.025 

Strain hardening exponent 

0.1 

0.095 

0.5 

Strain hardening modulus 

125.0 

1000.0 

2.0 

Melt temperature (degrees Kelvin) 

1,220 

1,710 

2,310 

Specific heat (Mbar-cm 3 per g-kilodeg Kelvin) 

0.00884 

0.00111 

0.00448 

Spall stress (Mbar) 

0.012 

0.028 

0.032 

Plastic failure strain 

1.0 

1.0 

1.0 


The first example problem involves the oblique (45 degree) impact of an 0.953 cm diam- 
eter aluminum sphere on a thin (0.1143 cm) plate, at 6.56 kilometers per second, and is 
representative of typical Whipple shield design problems. Figure 1 shows an element plot of 
the initial configuration of the projectile and target, while Figure 2 shows a particle plot of 
the simulation results at 6.6 microseconds after impact. The simulation results show good 
agreement with the experimental radiograph [17]. This problem was run on 16 processors of 
an IBM Regatta at two different mesh densities, with models composed of 0.67 million and 
3.20 million particles. The dimensions of the plate perforation were compared to determine 
the effect of model resolution on the simulation results. Table 3 shows that an increase in 
the particle count by nearly a factor of five produced only a small variation in the predicted 
dimensions of the plate perforation. 
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Fig. 3. Element plot of the initial 
configuration for the long rod impact 
problem. 


Fig. 4. Particle plot of the simulation 
results for the long rod impact problem. 


Table 3. Simulation results for the sphere impact problem 


Problem size 

Wall clock time 

Perforation width 

Perforation length 

(particles) 

(hours) 

(cm) 

(cm) 

0.673 million 

4.04 

1.90 

2.43 

3.196 million 

39.6 

1.81 

2.39 


The second example problem involves the oblique (73.5 degree) impact of an 0.767 cm 
diameter uranium alloy rod (L/D = 10) on a steel plate target (thickness 0.64 cm), at a 
rod velocity of 1.21 kilometers per second (target velocity was 0.217 km/s). This problem 
is representative of typical armor- antiarmor design applications. Figure 3 shows an element 
plot of the initial configuration of the projectile and target, while Figure 4 shows a particle 
plot of the simulation results at 100 microseconds after impact. The simulation results 
provided in Table 4 show good agreement with the corresponding experimental values [18] of 
residual rod length (5.55 cm) and residual rod velocity (1.07 km/sec). This problem was run 
on 16 processors of an IBM Regatta at two different mesh densities, with models composed 
of 1.56 million and 5.06 million particles. Table 4 shows that a factor of more than three 
increase in the particle count produced only a small variation in the simulation results. 


Table 4. Simulation results for the long rod impact problem 


Problem size 

Wall clock time 

Residual length 

Residual velocity 

(particles) 

(hours) 

(cm) 

(km/s) 

1.566 million 

13.8 

5.74 

1.09 

5.060 million 

74.4 

5.78 

1.09 


The third example problem involves the oblique (45 degree) impact of an aluminum shaped 
charge projectile (0.5062 cm diameter, 2.2046 cm length) on an aluminum plate target pro- 
tected by a dual plate aluminum debris shield. The projectile velocity was 11.0 kilometers 
per second, while the wall plate thickness (0.50 cm), shield thickness (0.25 cm), and total 
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Fig. 5. Element plot of the initial 
configuration for the dual plate shield 
problem. 


Fig. 6. Particle plot of the simulation 
results for the dual plate shield problem. 


standoff distance (12.0 cm) are representative of orbital debris shielding design applications. 
This problem has been designated a benchmark for use in the numerical analysis of spacecraft 
protection systems [19]. Figure 5 shows an element plot of the initial configuration of the sys- 
tem, while Figure 6 shows a particle plot of the simulation results at 150 microseconds after 
impact. Consistent with the corresponding experiment, the simulation results show bulging 
but no perforation of the wall plate. This problem was run on SGI Origin systems at two 
different mesh densities, with models composed of 4.27 million and 12.90 million particles. 
The smaller model required 56.8 wall clock hours on 256 (400MHz) processors, while the 
larger model required 332 wall clock hours on an average of 502 (600Mhz) processors. The 
two simulations indicated very similar wall plate damage. Figures 7 through 9 show details 
of the shield perforations and wall plate damage, in element plots made at 150 microseconds 
after impact. This problem illustrates that the hybrid numerical method used here is well 
suited to represent both the very general contact-impact dynamics illustrated in Figure 6 
and the large deformation plasticity illustrated in Figures 7 through 9. 

4. Parallel Speedup Tests 

The simulations described in the preceding section illustrate that three dimensional hy- 
pervelocity impact simulation can be very computer resource intensive. Hence the parallel 
performance characteristics of the relevant algorithm and numerical implementation is of 
considerable interest. This section discusses several algorithmic and implementation issues 
and presents measured parallel speedup data. Parallel implementation [20] of the hybrid 
method discussed here is based on the OpenMP standard [21], a set of portable compiler di- 
rectives which may be simply applied to achieve loop level parallelism. Although an OpenMP 
implementation does not manage distributed memory, the effects of nonuniform memory ac- 
cess on code performance can in some cases be observed and influenced, through changes in 
the initial data placement scheme. 

The hybrid nature of the present numerical formulation requires that both finite element 
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Fig. 7. Element plot of the simulation 
results for the dual plate shield problem. 


Fig. 8. Oblique view of the simulation 
results for the dual plate shield problem. 


based and particle based computation take place throughout the simulation. The Lagrangian 
finite element calculations involve an invariant neighbor set for each particle, known at the 
start of the calculation. This allows for a very efficient parallel implementation, and sug- 
gests the use of a material based domain decomposition technique on distributed memory 
systems. Unfortunately a Lagrangian description of contact-impact and fragmentation ef- 
fects, applied here and represented using the particle kinematics, involves a second (and 
time varying) neighbor set for each particle. Calculations involving this neighbor set are 
relatively inefficient, since a significant number of nominal neighbors identified by a search 
algorithm will turn out to be outside the effective mechanical or thermal interaction range. 
This portion of the calculation suggests the use of a geometry based domain decomposition 
technique, like that used in Eulerian codes, for distributed memory systems. However in 
this case, unlike Eulerian models, the membership of the contact-impact neighbor set for a 
given particle is time varying. This can present significant problems with load balancing and 
communications overhead. 

The aforementioned considerations are reflected in the results of performance testing on 
the code used here, which indicates that: (a) particle based calculations largely determine the 
required wall clock time, (b) most wall clock time is spent in the routines which calculate the 
densities and particle interaction forces, via summations over the contact-impact neighbor 
set, and (c) round-robin (or maximum entropy) initial data placement provides the best 
overall performance on nonuniform memory access systems. 

The authors provided in reference [3] speedup data for test cases run on as many as 128 
processors of an SGI Origin, applying a hybrid particle-element modeling methodology to 
problems as large as 0.5 million particles. As noted in the last section, more recent work has 
focused on much larger models, hence speedup tests have recently been performed at larger 
processor counts. These tests employed a 1024 processor SGI Origin and problem sizes as 
large as fifteen million particles. The test problem used for speedup measurements was the 
oblique sphere impact problem discussed in the section on validation simulations. Tables 5 
and 6 show the results of test cases run for two different model sizes (4.7 and 14.6 million 
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Fig. 9. Second oblique view of the simulation results for the dual plate shield problem. 


particles), at various processor counts, ranging from 256 to approximately one thousand. In 
both cases speedup is measured relative to the wall clock time required for the test problem 
on 256 processors, and efficiency is defined as the ratio of measured to ideal relative speedup. 
The wall clock times shown are those required to complete 100 time steps of the simulation. 

Table 5 shows good relative efficiency for the smaller test problem, when moving from 
256 to 504 cpus, but very poor performance (in fact an increase in wall clock time) when 
the cpu count is then increased to 1008. Note that in the 1008 processor simulation each 
cpu is allocated approximately 5,000 particles, a load level apparently too light to allow for 
efficient parallel performance. Table 6 provides results for the larger test problem, where 
a minimum load level of nearly 15,000 particles per processor is maintained. In this case 
parallel performance is significantly improved, with a relative efficiency of almost seventy 
percent measured as the processor count is increased from 256 to 976. The preceding results 
indicate that the numerical method and OpenMP implementation tested here can efficiently 
address rather large scale problems. A dependence of parallel efficiency on processor load, 
observed here, is not unusual for engineering applications [22], 


Table 5. Speedup measurements for a 4.7 million particle test problem 


Number of 
processors 

Particles 
per processor 

Wall clock time 
(hours) 

Speedup 

(relative) 

Efficiency 

(relative) 

256 

18,281 

0.7858 

1.000 

1.000 

504 

9,286 

0.4858 

1.618 

0.822 

1008 

4,643 

0.5008 

1.569 

0.400 
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Table 6. Speedup measurements for a 14.8 million particle test problem 


Number of 
processors 

Particles 
per processor 

Wall clock time 
(hours) 

Speedup 

(relative) 

Efficiency 

(relative) 

256 

57,031 

1.8439 

1.000 

1.000 

512 

28,516 

1.0391 

1.775 

0.887 

640 

22,813 

0.8803 

2.095 

0.838 

976 

14,959 

0.7031 

2.623 

0.688 


5. Conclusion 

A number of new numerical methods based entirely or in part on particle kinematics 
are currently under development. They emphasize the fact that some important engineer- 
ing problems are dominated by noncontinuum physics, such as fragmentation and contact- 
impact. Such problems may not fit easily into a classical continuum mechanics modeling 
framework. The present paper has described recent work aimed at extending and validating 
a hybrid numerical method for hypervelocity impact simulation. The results suggest that the 
method is accurate, numerically robust, and suitable for large scale parallel computation. 
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Ravishankar Shivarama 3 and Eric P. Fahrenthold 4 
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A combination of Euler parameter kinematics and Hamiltonian mechanics provides a rigid 
body dynamics model well suited for use in strongly nonlinear problems involving arbitrarily 
large rotations. The model is unconstrained, free of singularities, includes a general poten- 
tial energy function and a minimum set of momentum variables, and takes an explicit state 
space form convenient for numerical implementation. The general formulation may be spe- 
cialized to address particular applications, as illustrated in several three dimensional example 
problems. 


INTRODUCTION 

A variety of rigid body dynamics modeling problems demand consideration of very large 
rotations. Some of the best known examples involve aircraft [1] and spacecraft [2,16], al- 
though the analysis of large rotation dynamics is of generic interest in a wide range of 
applications, including mechanism and machine theory [3,4] and molecular dynamics [12]. 
Most models of rigid body dynamics problems employ Euler angles [5]. Such formulations 
lead to equations of motion which are unconstrained, but which contain singularities [1]. 
Other singular three parameter methods have been developed, including for example those 
of Laning-Bortz-Stuelpnagel [11] and Rodriguez [13]. The presence of singularities in all these 
methods has motivated the development of alternative four parameter modeling schemes [11], 
including Euler parameters [15]. Such formulations replace the three Euler angles with four 

3 Graduate research assistant 

4 Professor 
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parameters and an algebraic constraint. This avoids the Euler angle singularities but leads 
nominally to a system level model in differential-algebraic form. 

In an attempt to avoid both singular equations of motion and differential-algebraic sys- 
tems, several authors have presented reformulations of Euler parameter based models, for 
use in three dimensional rigid body dynamics problems. Chang et al. [2], Nikravesh and 
co-workers [7,8,9,10], and Vadali [17] present alternative formulations based on Lagrange’s 
equations. Since Lagrange’s method defines the solution as a path in configuration space, 
and since the Euler parameters are taken as generalized coordinates, this approach starts 
with a differential system of order eight (four second order equations for the rotational dy- 
namics of a single rigid body) augmented with a single algebraic constraint. Nikravesh and 
co-workers begin from this starting point and proceed to find a closed form solution for the 
Lagrange multiplier associated with the algebraic constraint, resulting in an unconstrained 
formulation of order eight. They do not include a potential energy function in the system 
Lagrangian. Similar results are obtained by Vadali. Proceeding in a different manner, Chang 
et al. introduce as quasi-velocity variables the rigid body angular velocities in the body fixed 
frame, and project the original order eight Lagrange equations onto an order seven subspace. 
In the process they eliminate the unknown Lagrange multiplier. 

As an alternative to Lagrange’s equations, a Hamiltonian formulation of rigid body dy- 
namics with Euler parameters has been proposed by Morton [6]. However his final formula- 
tion is of order eight, and includes a superfluous momentum variable as well as a ‘generally 
arbitrary’ unspecified scalar parameter. It appears that no previous work has attempted to 
revise or improve upon the Morton formulation. 

The usefulness of formulations based on Hamilton’s canonical equations is well recognized 
[3]. They offer an explicit state space description of system dynamics problems which is: 
(a) convenient for numerical integration, (b) well suited for coupling to automatic control 
system models, and (c) energy based and hence providing clear physical insight. Recognizing 
these strengths, a revision and extension of existing Hamiltonian formulations for rigid body 
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dynamics is of generic interest. The present paper presents such work, deriving unconstrained 
Hamilton’s equations for the three dimensional dynamics of a rigid body in terms of Euler 
parameters, and hence suitable for use in simulations involving arbitrary rotational motion. 
The derivation avoids any requirement to determine the Lagrange multiplier associated with 
the Euler parameter constraint. No arbitrary parameters are introduced, and the final 
rotational formulation is of order seven. A general potential energy function and nonpotential 
virtual work effects are included in the model. Validation and application of the method is 
illustrated here in several three dimensional example problems. 

KINEMATICS 

This section defines the kinematic variables of interest, and recalls a number of well known 
kinematic relations [1], for use in succeeding sections. 

The position and orientation of an arbitrary rigid body is described here in terms of seven 
generalized coordinates, namely the Cartesian components of the center of mass vector ( c ) 
and a four component vector of Euler parameters ( e ) 

c = [ ci c 2 c 3 ] T , e = [ e 0 ei e 2 e 3 } T (1) 

Knowledge of the Euler parameters determines a (nonunique) set of Euler angles (0, 9, ij>) 
for the body, associated with a 3-1-3 rotation sequence, via the relations 


0 

= tan *( — ) + tan x ( 


(2) 


= fa77, _1 ( — ) — tan~ l { 

A 

(3) 


e 0 

ei 


9 

= 2 sin~ 1 ( \/ ei 2 + e 2 2 

) 

(4) 


The Euler parameters define an orthogonal rotation matrix (R) which relates the vector of 
components (a) of a first order tensor described in a fixed Cartesian coordinate system to a 
corresponding vector of components (a) described in a body fixed co-rotating frame, using 

a = R a (5) 
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where 


with 



R 

= EG 

T 



— ei 

e o ~ 

e 3 

e 2 

E = 

-e 2 

e 3 

e 0 

-ei 


. -e 3 

e 2 

ei 

e 0 . 



-ei 

e 0 

e 3 

-e 2 

G = 

-e 2 

— e 3 

e 0 

ei 


-e 3 

-e 2 

-ei 

e 0 


( 6 ) 

( 7 ) 

( 8 ) 


The four Euler parameters are not independent, and satisfy the constraint equation 

e T e = 1 (9) 

which then implies 

G G r = I (10) 

where I is an order three identity matrix. In addition, G and e and their time derivatives 
satisfy the identities 

G e = 0, G e = — G e (11) 

The kinematic equations [1] which relate the time derivatives of the Euler parameters 
to the components of the angular velocity vector ( uj ) of the rigid body, expressed in the 
body-fixed co-rotating frame, are 

w = 2Ge, e = ^ ^ T ^ (12) 

Finally note that the skew-symmetric matrix 12, with axial vector ay which satisfies 

f2v = wxv (13) 

for any vector v, is related to the Euler parameters and their time derivatives by 

n = 2 G G t = — 2GG t (14) 
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The next section defines kinetic and potential energy functions and hence the Hamiltonian 


for the system of interest. 


KINETIC AND POTENTIAL ENERGY 

The complementary kinetic energy for the rigid body may be expressed as 

T* = - m c T c + - uj 1 J lo 
2 2 


(15) 


where m is the mass and J is a constant matrix of components for the moment of inertia 
tensor of the body, referred to the co-rotating frame. In terms of the Euler parameters and 
their time derivatives, 

r = ^mc T c + 2e T G T JGe (16) 

which has the form T* = T* ( c, e, e ). It follows that the generalized momenta are 


dT* 

P = ~ rp = rn c , 

a c 


dT* T 

g= — =4G r JGe 
ae 


(17) 


Note that the identities (10) through (12) require 


g = 2 G t J uj , 


uj = - J^Gg 
2 B 


(18) 


Since the complimentary rotational kinetic energy may also be expressed, using equation 
(11), as 

r; ot = 2e r G r JGe (19) 

then the Euler parameter dependence of T* defines the partial derivative 


dT* 

k = =4G T JGe 

de 

The kinetic energy of the body is defined via the Legendre transform 


( 20 ) 


T = p J c + g T e - T 


so that the preceding results lead to the canonical form T = T ( p, g, e ) which is 

T- ( p'p I ‘ g''G r .l r G g 


( 21 ) 


( 22 ) 
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and require that (see the appendix) 


, dT 


( 23 ) 


For the mechanical systems considered here, the potential energy function has the general 
form 

V — V(c, e) (24) 

and the system Hamiltonian is 

H = T + V (25) 

The next section introduces a virtual work expression, to account for nonpotential effects. 


NONPOTENTIAL VIRTUAL WORK 

The quasi-coordinates q associated with the co-rotating components of the angular veloc- 
ity vector are defined by 

q = u) (26) 

In terms of the latter coordinates, and the center of mass coordinates, the nonpotential 

virtual work due to the imposed forces f (t) and torques T(i) is 

5W nc = f(t) T hc + T(f) T hq (27) 

Since 

hq = 2 G he (28) 

the virtual work expression which defines the generalized nonpotential forces is 

8W nc = f (t) T hc + 2 [G t T (t) ] T he (29) 

Note that damping effects may contribute additional terms to the nonpotential virtual work, 
in which case the nonpotential forces may depend on the generalized coordinates and ve- 
locities. In addition the presence of nonholonomic constraints may introduce terms which 
depend on unknown Langrange multipliers. The last problem discussed in the examples 
section illustrates the effects of both damping and nonholonomic constraints. 
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The next section derives Hamilton’s equations for the system. 


HAMILTON’S EQUATIONS 

The system Hamiltonian has the form H = H ( p, c, g, e ) and the canonical Hamilton’s 
equations are 

dH dH 

P = -^ + f p , c = — (30) 


and 


g = 


dc 

dH 


f 9 . 


dp 

dH 


e = 


(31) 


de ’ d g 

where f p and f 9 are the nonpotential generalized forces associated with the virtual work and 
any applied constraints. The Euler parameter constraint has the rate form 


e 1 e = 0 


(32) 


Introducing a Lagrange multiplier A, the latter constraint combines with the virtual work 
expression to yield 


f p = f(t) 

f 9 = 2 G T T(f) + Ae 

so that for the derived Hamiltonian the momentum balance equations are 

dV 

p = -^ + f(<) 

dV 

g = + k + 2G T T(t) + Ae 

ae 


(33) 

(34) 


(35) 

(36) 


The last equation includes an unknown Lagrange multiplier and a superfluous momentum 
variable. These variables are eliminated by introducing the three-momentum vector 


h = Jw 


(37) 


whose time derivative is 


h= ^Gg + ^Gg 


(38) 
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With equations (11) and (36) this yields 



( 39 ) 


or 


11 1 dV 

h = --n h + - Gk + T (f)--c;^ 


(40) 


which eliminates both A and the four component momentum vector. 


The final unconstrained Hamiltonian model is 


P 


h 



(41) 


(42) 


c 



(43) 


e 




(44) 


Given a potential function and the virtual work, the preceding explicit equations may be 
integrated to simulate the system response. 

As outlined in the introduction, the rigid body dynamics formulation derived here com- 


nian and Lagrangian methods generally simplify the model formulation process, in particular 
when geometric nonlinearities are important. Such nonlinearities arise for example when 


normally most convenient for numerical integration. Euler parameters offer a singularity 
free description of rotational displacements, and are therefore preferred over Euler angle 
models in large rotation applications. The cost is of course the need to integrate an addi- 
tional state equation for each rigid body, since the Euler parameters are a quaternion. The 
present combination of Hamiltonian mechanics and Euler parameter kinematics is therefore 
of most interest in the formulation and numerical integration of models for strongly nonlin- 
ear mechanical systems. The model developed here is unique in its combination of features: 


bines the advantages of Hamiltonian mechanics and Euler parameter kinematics. Hamilto- 


large rotations or hyperelastic devices are of interest. As compared to Lagrangian meth- 
ods, Hamiltonian methods offer an explicit state space description of the system dynamics, 
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unconstrained, free of singularities, incorporating a general potential energy function, em- 
ploying a minimum set of momentum variables, and taking an explicit state space form 
convenient for numerical implementation. 

EXAMPLE PROBLEMS 

Application of the Hamiltonian formulation developed here is illustrated in four examples. 
The first example compares an Euler angle based model of a rotating disk problem to the 
present Hamiltonian formulation, both to validate the present approach and to illustrate 
in a simple case an Euler parameter description of rotational displacement. The second 
example solves a classical rigid body dynamics problem for the three dimensional motion of 
a torque-free body, for comparison to the published numerical solution of Morton [6] and to 
the partial analytical solution of Thompson [16]. This problem again validates the present 
approach, and compares a discontinuous Euler angle description of three dimensional rigid 
body motion to the continuous Euler parameter characterization adopted here. The third 
example models a spinning top in a uniform gravitational field, a problem described in many 
advanced dynamics tests, and validates the present formulation in a three dimensional rigid 
body motion involving both translation and rotation. Here a numerical solution for the last 
cited problem is obtained using a time step identical to that employed by Simo and Wong 
[14], but without resort to their symplectic integration algorithm. The last example considers 
a problem of practical importance in the design of gyroscopic seekers, namely the motion of 
a freely precessing body with a viscous ring nutation damper [2], This example calls for the 
application of nonholonomic constraints. Here we develop an explicit state space model, as 
compared to the implicit Lagrangian formulation [2] of Chang et al. 

The first example models the free vibration of a rigid circular disk of radius r, rotating 
about a fixed point, and attached to a linear spring of stiffness k (see Figure 1). The potential 
energy function is 

V = ^ ky 2 p = l,kr 2 [ 2eie 2 + 2e 0 e 3 ] 2 (45) 
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where y v denotes the vertical displacement of the point of attachment of the spring, measured 
in a Cartesian coordinate system whose origin lies at the center of the disk. Note that the 
indicated Euler parameter dependence of the potential energy is obtained using equation 
(5). Assuming the model parameters and initial conditions listed in Table 1, the motion was 
simulated by integration of a Newtonian model based on the Euler angle cj) 

J(j)+^kr sin{ 20) = 0 (46) 

and by integration of the Hamiltonian relations (42) and (44). The two computed results for 
the angular momentum are compared in Figure 2, showing excellent agreement of the Euler 
angle based and Euler parameter based solutions. Figure 3 shows the computed variation of 
the Euler parameters with time. 

The second example models the torque free motion of a rigid body, for the inertial prop- 
erties and initial conditions listed in Table 2. A partial analytical solution for this classic 
problem is known and can be expressed in terms of elliptic functions [6,16]. Figures 4 and 
5 shown the time variation of the angular momenta and Euler parameters computed using 
the present Hamiltonian formulation. Figure 6 plots the implied Euler angles, emphasizing 
the discontinuous nature of the latter variables. Table 3 shows excellent agreement of the 
analytical and numerical solutions for the amplitudes and periods of the angular momenta, 
in three dimensional motion. 

The third example models the translational and rotational motion of a spinning top in 
a uniform gravitational held. This problem is described in many advanced dynamics texts 
[4], and is used by Simo and Wong [14] to evaluate their symplectic numerical integration 
scheme. The potential energy for the system is 

V = W z c (47) 

where W is the weight of the top and z c is the vertical coordinate of the center of mass. The 
simulation parameters and initial conditions for the problem are listed in Table 4. Figures 7 
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and 8 show numerical results for the normalized angular momentum components and center 
of mass coordinates, obtained by integration of Hamilton’s equations (41) through (44), 
using a fourth order Runge-Kutta method. Table 5 compares an approximate analytical 
estimate of the nutation and precession frequencies for the top, provided by Goldstein [4], 
to the present numerical results. The present numerical results are identical to those plotted 
by Simo and Wong [14], and are obtained using the same time step, but without resort to 
their symplectic integration scheme. 

The fourth example considers the rotational motion of a rigid rotor damped by a partially 
filled mercury ring damper. The reader is referred to Chang et al. [2] for a detailed discussion 
of this problem, and its application in the analysis of gyroscopic seekers. We focus here on 
the formulation of a dynamic model for the system analyzed in reference [2], The paragraphs 
which follow develop an explicit Hamiltonian model for this system, an alternative to the 
implicit Lagrangian model of Chang et ah, adopting their stipulated assumptions on stored 
energy functions, energy dissipation, and kinematic constraints. 

The rotor is modeled as a rigid circular cylinder with a fixed center of mass located at 
the origin of a global XYZ coordinate system. The partially filled mercury ring damper 
is a cylinder of mean radius R , co-axial with and external to the rotor, with a centroid 
displaced a distance L along the rotor axis from the rotor center of mass location. Body- 
fixed coordinate systems for the rotor (xyz) and ring (uvz) are co-located at the centroid of 
the damper, where the z direction is aligned with the rotor axis. The partial mercury ring is 
free to rotate about the rotor axis, subject to a damping torque which is linear in the axial 
angular velocity difference between the rotor and the ring, but is otherwise constrained to 
move with the rotor. Hence the orientations of the body-fixed axes systems which co-rotate 
with the rotor and the ring differ only by an angle (3, which describes the axial rotation of 
the ring with respect to the rotor. 
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The assumed complimentary kinetic energy for the system is [2] 


r = \ ^ T J T « + \ (48) 

where U3 and are angular velocities for the rotor and ring and 


■ •/, 
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1 

CO 


are constant moment of inertia matrices for the rotor and ring. All four quantities are 
described in the respective rotor and ring body-fixed co-rotating coordinate systems. 

The assumed potential energy for the system, due to the gravitational potential of the 
mercury, is [2] 

V = -m g c R T B T r c , r c = [R sin{^/2)/{^/2), 0, L} T (50) 


where m is the mass of the mercury, g c is a constant gravity acceleration vector described 
in the fixed XYZ frame, R is the rotation matrix of equation (6), whose Euler parameters 
(e) refer to the rotor-fixed frame, r c is a constant vector which locates the mercury center 
of mass, 7 is the angle which subtends the mercury arc (symmetric about the u axis), and 
B is an orthogonal matrix which defines the transformation of vector components from the 
rotor-fixed to the ring-fixed frame 


B 


cos(/3 ) sin(fl) 0 
— sin(/3 ) cos(fi) 0 

0 0 1 


(51) 


Note that V = V(e,/3). The virtual work for the system, due to damping at the ring-rotor 
interface, is [2] 

5W = -C d R 2 ft 6P (52) 


where Cd is an empirical dimensionless damping coefficient. 

Given the preceding modeling assumptions, Hamilton’s equations for the rotor and ring 
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system are 


h = 

-fih-^G^ + T 

2 oe 

(53) 

h m 

+ T m 

(54) 

e = 


(55) 

0 = 

dV m 
~~dp + p 

(56) 


where h and h m are angular momenta for the rotor and ring 

h — J uj , h m = 3 m uj m (5 7) 

and T. T m , and Tp are nonpotential forces due to damping and kinematic constraints. Note 
that the degenerate form of Hamilton’s equation for (3 is due to the fact that the latter 
generalized coordinate, which appears in the potential energy function, is not associated 
with a corresponding generalized momentum variable. The kinematic constraints are [2] 

Umi = cos{(3) uii + sin(P) lu 2 , Cc ; m2 = — sin((3 ) LUi + cos{(3 ) u 2 , (58) 

and 

P = U>m3 - ^3 (59) 

They quantify the aforementioned modeling assumption that the ring moves relative to the 
rotor only in axial rotation. 

An explicit state space model may be obtained by application of the constraints, as follows. 
Introducing a Lagrange multiplier /i for the constraint (59), and accounting for the virtual 
work, requires 

T /3 = n - C d R 2 (3, T = n c, T m = -fx c (60) 

where c denotes the vector [0,0, 1] T . The degenerate Hamilton’s equation for (3 therefore 
determines the Lagrange multiplier as 

8V 

H= C d R 2 (3+ — (61) 
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Hamilton’s equation (54) for the ring angular momentum may now be written in the form 


<jJ r 


— J™ ^m.Jr 




h Jr 


(62) 


Since the constraints specify both (3 and the first two components of the ring angular velocity 
vector as functions of the set (f3,coi,U2,^m 3 ), the third of equations (62) is an evolution 
relation for the unknown uJm 3 . Combining the third of equations (62) with the constraint 
equation (59), the constitutive relations (57), and Hamilton’s equations (53) and (55), the 
result is an explicit state space model of order nine for the ring-rotor system. The final state 
equations are 


h = 


e = 


P = 


^m3 — 


_f!h_ 2 G ^ + 


C, U 1 


h 3 \ dV 

“ ” 3 ~TJ + W. 


-C T J- lh 

h'i 

W m 3 V~ 

"3 

^ (j m 


C d R 2 


h 3 \ dV' 


c T Jjc 


where 


f2 — fi(h), G — G(e), Q, m — Q m (u: m ) 


(63) 

(64) 

(65) 

( 66 ) 

(67) 


with 


cVmi = cos({3) ^ + sin(P) 

J 1 J 2 


w m2 = - sin(/3 ) ^ + cos(f3) (68) 

J 1 J 2 


Note that implicit model of reference [2] is also of order nine but employs a different set of 
state variables. 


CONCLUSION 

The present paper has derived and applied a new Hamiltonian formulation of rigid body 
dynamics problems, based on Euler parameter kinematics. Euler parameter kinematics pro- 
vide a singularity free description of three dimensional rigid body motion, accommodating 
arbitrarily large rotations. When combined with Hamiltonian mechanics, the result is an 
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energy based modeling approach well suited to address problems with complex geometric 
nonlinearities. As compared to previous work, the formulation derived here offers a unique 
combination of features. It avoids the introduction of algebraic constraints and unspecified 
parameters, includes a general potential energy function, incorporates a minimum set of 
momentum variables, and takes an explicit state space form convenient for use in control 
related applications. 
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APPENDIX 


A complimentary kinetic energy expression (T*) with the functional form 


T* = T*(e,f), f = e 


( 69 ) 


has the total differential 


f)T* T 

dT* = g T di + — de , 
de 


g = 


dT* 

w 


(70) 


The corresponding kinetic energy function (T) is determined by the Legendre transform 


T = g T f - T* 


(71) 


and has a total differential defined by 


BT* t 

dT = g T di + f r dg — dT* = f r dg — de (72) 

oe 


as well as the canonical form 


rrTi 3 T t , 3T t , 
dT = — dg + — de 
a g oe 


It follows that 


dT dT _ dT* 

dg ’ de de 


(73) 


(74) 
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Parameter 

value 

Mass moment of inertia (kg m 2 ) 

J = 2 

Radius of the disk (m) 

r — 1 

Stiffness of the spring ( N/m ) 

k = 10 

Initial displacement ( degrees ) 

0 = 30 

Initial momentum (kg m 2 rad/s) 

h = 0 


Table 1: Model parameters and initial conditions for the first example problem 


Parameter 

value 

Mass moments of inertia (kg m 2 ) 
Initial Euler parameters 
Initial momenta (kg m 2 rad/ s) 

J\ = 400, J 2 = 307.808385, J 3 = 200 
e 0 = 1, ei = e 2 = e 3 = 0 
h i = 346.4101616, h 2 = 0, h 3 = -200 


Table 2: Model parameters and initial conditions for the second example problem 


Variable 

exact solution 

numerical solution 

Magnitude of h\ (kg m 2 rad/s) 

346.4102 

346.38 

Magnitude of h 2 (kg m 2 rad/ s) 

365.447 

365.44 

Magnitude of h 3 (kg m 2 rad/s) 

200.0 

199.975 

Period of hi (s) 

9.3393 

9.35 

Period of h 2 (s) 

18.6786 

18.68 

Period of h 2 (s) 

18.6786 

18.68 

Minimum of h 3 (kg m 2 rad/ s) 

162.6296 

162.6342 


Table 3: Exact versus numerical results for the second example problem 


61 


Parameter 

value 

Weight ( kgm/s 2 ) 

Mass moment of inertia (kg m 2 ) 
Initial Euler parameters 
Initial angular momenta (kg m 2 rad/s) 

W = 20 

Ji = 5, J 2 = 5, J 3 = 1 

e 0 = cos(0.15), e\ = sin(0.15), e 3 = 0, e 4 = 0 
hi = 0, h 2 = 0, h 3 = 50 


Table 4: Simulation parameters and initial conditions for the third example problem 


Variable 

analytical approximation 

numerical simulation 

Nutation frequency (rad/s) 

10.00 

9.24 

Precession frequency (rad/ s) 

0.40 

0.4136 


Table 5: Approximate analytical versus numerical results for the third example problem 
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Figure 1. First example problem, rotating disk with a translational spring suspension 

Figure 2. First example problem, comparison of Euler angle based and Hamiltonian solu- 
tions for the angular momenta versus time; the computed Hamiltonian component h z agrees 
with the Euler angle solution (diamond symbols) obtained using equation (46), while the 
Hamiltonian components h x and h y are identically zero 

Figure 3. First example problem, numerical solution for the Euler parameters versus time; 
note that the Euler parameters e± and e2 are identically zero, and along with the computed 
solutions for eo and e 3 determine the nonzero Euler angle in accordance with equation (2) 

Figure 4. Second example problem, torque free motion of a rigid body, numerical solution 
for the angular momenta versus time 

Figure 5. Second example problem, torque free motion of a rigid body, numerical solution for 
the Euler parameters versus time; note that the Euler parameters are continuous functions 

Figure 6. Second example problem, torque free motion of a rigid body, computed Euler 
angles versus time; note the discontinuities in the Euler angles 

Figure 7. Third example problem, translation and rotation of a spinning top, numerical 
solution for the normalized components of the angular momentum versus time 

Figure 8. Third example problem, translation and rotation of a spinning top, numerical 
solution for the center of mass position versus time 
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Figure 1: First example problem, rotating disk with a translational spring suspension 



time in secs 

Figure 2: First example problem, comparison of Euler angle based and Hamiltonian solu- 
tions for the angular momenta versus time; the computed Hamiltonian component h z agrees 
with the Euler angle solution (diamond symbols) obtained using equation (46), while the 
Hamiltonian components h x and h y are identically zero 
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time in secs 


Figure 3: First example problem, numerical solution for the Euler parameters versus time; 
note that the Euler parameters e± and e 2 are identically zero, and along with the computed 
solutions for eo and e 3 determine the nonzero Euler angle in accordance with equation (2) 
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400 



Figure 4: Second example problem, torque free motion of a rigid body, numerical solution 
for the angular momenta versus time 
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Figure 5: Second example problem, torque free motion of a rigid body, numerical solution for 
the Euler parameters versus time; note that the Euler parameters are continuous functions 
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Figure 6: Second example problem, torque free motion of a rigid body, computed Euler 
angles versus time; note the discontinuities in the Euler angles 
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h/llhll 



Figure 7: Third example problem, translation and rotation of a spinning top, numerical 
solution for the normalized components of the angular momentum versus time 
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Center of mass location 



Figure 8: Third example problem, translation and rotation of a spinning top, numerical 
solution for the center of mass position versus time 
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A KERNEL FREE PARTICLE-FINITE ELEMENT 
METHOD FOR HYPERVELOCITY IMPACT SIMULATION 


Young-Keun Park 3 and Eric P. Fahrenthold 4 

Department of Mechanical Engineering, 1 University Station C2200 
University of Texas, Austin, TX 78712, USA 

An improved hybrid particle-finite element method has been developed for the simulation of 
hypervelocity impact problems. Unlike alternative methods, the revised formulation computes 
the density without reference to any kernel or interpolation functions, for either the density 
or the rate of dilatation. This simplifies the state space model and leads to a significant 
reduction in computational cost. The improved method introduces internal energy variables 
as generalized coordinates in a new formulation of the thermomechanical Lagrange equations. 
Example problems show good agreement with exact solutions in one dimension and good 
agreement with experimental data in a three dimensional simulation. 

KEYWORDS: particle methods, finite element methods, impact simulation 

INTRODUCTION 

Studies of hypervelocity impact phenomena are motivated by a variety of science and en- 
gineering applications [1]. Examples include scientific research on planetary impacts [2] and 
equations of state [3] and engineering research on the design of spacecraft shielding [4] and 
kinetic energy penetrators [5] . The proceedings of a recent international symposium [1] show 
that the use of computer simulation in this field is increasing, as improvements in numerical 
methods and computing power make it possible to address problems of greater complexity 
and larger scale. Simulation is of particular importance, as an adjunct to experimental work, 
when material costs are high [6] or when impact velocities beyond the range of light gas guns 
are of interest [7]. 

Simulation work in this field has applied a number of different numerical methods, based 
on continuum mechanics, particle dynamics, or mixed kinematic schemes. Continuum meth- 
ods [8] employ either an Eulerian hydrodynamic [9,10] or a Lagrangian finite element [11] 

3 Graduate research assistant 

4 Professor, corresponding author, phone: (512) 471-3064, email: epfahren@mail.utexas.edu 
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approach, or some Arbitrary Lagrangian-Eulerian (ALE) based generalization of these tech- 
niques [12,13]. A large majority of particle codes employ a smooth particle hydrodynamics 
(SPH) technique [14,15,16], although some alternative particle based methods have been 
proposed [17]. Some disadvantages of pure continuum or pure particle based methods [18] 
have motivated the development of mixed continuum-particle formulations [4,19,20]. The 
most widely used mixed method is a coupled particle-finite element technique [19]. This 
technique initializes distinct material regions with either SPH particles or or Lagrangian 
finite elements, then quantifies subsequent material interactions using a particle-to-surface 
contact-impact algorithm. The alternative coupled particle-element method of Johnson and 
co-workers [20] maps damaged or failed elements into particles, and again quantifies particle- 
element interaction using a special contact algorithm. Both these methods are subject to 
tensile instability and numerical fracture problems. 

The alternative mixed method of Fahrenthold and co-workers [21] is based on a hybrid 
particle-finite element formulation. This method in not subject to tensile instability and 
numerical fracture problems and eliminates the requirement for special treatment of particle- 
to-element contact-impact. It avoids both the mass diffusion problems of Eulerian methods 
and the mass and energy discard associated with Lagrangian element erosion algorithms. It is 
labeled a hybrid (versus coupled) method since it introduces both elements and particles for 
all material control volumes, then employs the elements and particles in tandem to represent 
distinct physics. The particles model all inertia, contact-impact, and thermomechanical 
response in compressed states, while the elements model tension and elastic-plastic shear. 
The method incorporates both ellipsoidal particles and time varying particle volumes and as a 
result can represent large density variations with relatively small neighbor counts. Previous 
work employing nonspherical evolving kernels has been rather limited and most particle 
simulations represent high densities using spherical particles, a fixed contact length, and 
relatively large neighbor sets. 

The preceding formulation combines a true Lagrangian description of material strength 
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effects with a general particle based model of contact-impact dynamics, and has been vali- 
dated in simulations of impact experiments conducted at velocities ranging from one to ten 
kilometers per second [22], In the hypervelocity impact regime, where large strain plastic- 
ity, perforation, fragmentation, melting, and complex multi-structure contact-impact effects 
are often present, this formulation provides a particular combination of advantageous fea- 
tures not offered by alternative numerical methods. The present paper describes an im- 
proved hybrid particle-finite element method, modifying the formulation of Ravishankar 
and Fahrenthold [21] in two respects. First, the density is determined by integrating non- 
holonomic constraints imposed on the system level thermomechanical model. Second, the 
entropy states used previously to model the thermal domain are replaced by particle inter- 
nal energies. These modifications simplify the method, reduce its computational cost, and 
incorporate equations of state expressed in standard functional form. 

Unlike alternative methods, the revised formulation eliminates entirely the use of kernel 
or interpolation functions to represent the density or rate of dilatation fields. The density 
evolution equations are developed by direct reference to large deformation kinematics, avoid- 
ing any requirement to specify the functional dependence of an interaction potential on the 
particle coordinates. The latter task has proven to be quite difficult in an SPH context and is 
a principal focus of the general particle dynamics literature. The revised method introduces 
the use of internal energy variables as generalized coordinates in a new thermomechanical 
formulation of the discrete Lagrange equations. This avoids the requirement to construct 
Legendre transforms of the internal energy function, in order to express the dependence of 
pressure and temperature on entropy. 

The present paper is organized as follows. First the particle and element kinematics are 
defined, followed by the kinetic co-energy and thermomechanical potential energy functions 
for the particle-element system. Second the evolution equations for the density are developed, 
followed by the evolution equations for the plastic and damage variables, all of these relations 
representing nonholonomic constraints on the system level model. Third the numerical 
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viscosity and numerical heat diffusion models are introduced, and evolution equations for 
the internal energy state variables are described, the latter states serving as generalized 
coordinates in a thermomechanical Lagrangian formulation. Fourth the discrete Lagrange 
equations for the particle-element system are derived, taking an explicit state space form 
convenient for numerical implementation. Finally application of the method is illustrated 
in one dimensional problems with exact solutions and in three dimensional simulations of 
representative hypervelocity impact problems. 


PARTICLE KINEMATICS 

The inertia of the modeled system is represented by a collection of n ellipsoidal particles, 
with the mass of particle i and h^\ h^\ the half-lengths of its major axes. The 
position and orientation of each particle is determined by its center of mass position vector 
( cW ) and an Euler parameter [23,24] vector ( ) 


c 


(0 


r „(*) „(*) Ji) 1 T 

lb 4 b J , 


e« = 


,W Ji) JO JO 1 T 


( 1 ) 


where a superscript T denotes the transpose. 

It is convenient to note here certain properties of Euler parameters, and to cite a number 
of well known [23] kinematic relations associated with their use. The Euler parameters 
provide a singularity free description of arbitrary particle rotations. They define a rotation 
matrix (R^) for each particle 
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which relates vector components v described in a fixed global Cartesian coordinate system 
to corresponding components v described in a co-rotating system aligned with the particle 
major axes, using 

v = R w v (5) 

The Euler parameters and their time derivatives are related to the angular velocity vector 
of the particle ( c*/ 1 ' ), described in the co-rotating frame, by 

e w = ^ G wr u> (i) (6) 

Similarly the antisymmetric matrix with axial vector which satisfies 

O w v = w (i) x v (7) 

for all vectors v, is related to the Euler parameters and their time derivatives by the relations 

n (i) = 2 G (i) G (i)T = -2 G (i) G (i)r = R (l)T R (i) (8) 


For ellipsoidal particles it is convenient to describe the separation distance of the mass 
centers for particles i and j using the ellipsoidal coordinate 


C (iJ) = (c (i) - c 


0)) T Hd) ( c (0 _ c^) 


(9) 


defined in the co-rotating system of particle j using 

H W = Rb>H^')RW T 


( 10 ) 




2 (3h[ j) 0 0 

0 2 f3h { 2 j) 0 

0 0 2 fihf 


( 11 ) 


where the constant [3 allows for close packing at the reference density. The time derivative 
of this ellipsoidal coordinate , defined for i ^ j, is 


£(m) _ 




jjO') r (*j) J f .fid) (H^)r(* J ) x f-h'dy 


u> 


(j) 


( 12 ) 
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where 


r (*j) _ c (*) _ c 0) 




(13) 


and may be used to quantify the rate of compression for an array of ellipsoidal particles. 

The preceding results will be used in later sections to account for rotational inertia and 
particle kinematics not present in the vast majority of particle models, which assume a 
spherical particle geometry. 

FINITE ELEMENT KINEMATICS 

This section describes the finite element kinematics employed in the present paper. The 
elements used here are eight noded hexahedra, well known and described in detail by Hal- 
lquist [11] and others. Since all inertia effects are represented by the particles, no mass 
matrix is defined. 

Each structure in the model is subdivided into uniform hexahedra with orthogonal faces, 
with ellipsoidal particles located at each node and at the centroid of each element. The 
center of mass coordinates for particles located at element vertices are also nodal coordi- 
nates for the hexahedra, and are used to compute the shearing strain. The center of mass 
coordinates for particles located at the element centroids are used, in combination with the 
nodal coordinates, to define six subelements for each hexahedron. The volumes of these 
subelements are used to compute interparticle tension forces. 

The following Lagrangian finite strain deformation measures [25] are used in the stored 
energy functions for the elements, associated with tension and shear, and in the plastic 
constitutive relations. The shear strain for element j is 



(14) 


where 



and F-T is an element deformation gradient computed using one point integration [11], The 
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clastic shear strain is defined as 


E e 0') — — E p(j ') 


(16) 


where (E p -^) is a plastic stain tensor whose flow rule satisfies the isochoric plastic deformation 
constraint 


tr = 0, K" ,J I = i - 1) (17) 


( 17 ) 


The subelement Jacobians are denoted by J^ ,k \ where the index k designates one of six 


subelements for the jth hexahedron. 

KINETIC CO-ENERGY AND POTENTIAL ENERGY 

An energy method (Lagrange’s equations) is adopted here, to facilitate the systematic 
integration of diverse particle and element based modeling concepts. The stored energy 
functions considered here are a kinetic co-energy function for the particles, an internal energy 
function for the particles, and element potential energy functions which account for tension 
and shear. Damage variables are introduced to model element failure in a thermodynamically 
consistent fashion. Constitutive assumptions different from those adopted here may be 
introduced without change to the underlying methodology. 

The system kinetic co-energy is the sum of the particle co-energies 


n 



(18) 


where T*W is the co-energy for particle i, due to translation and rotation 



2 2 


(19) 


with Jh) a constant moment of inertia matrix described in the co-rotating particle frame. 
The system kinetic co-energy function defines the generalized momenta 




( 20 ) 


where pW and are translational and angular momentum vectors for the Ah particle. 
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The system potential energy has the general form 


1/ = u w + y v o U) ^ U) + Y v o (j,k) ^ {j,k) 

i = 1 j = 1 j = 1 fc = 1 


( 21 ) 


where is the total internal energy for particle i and the pressure (PW) and temperature 
($W) are described by an equation of state [26] with functional form 

pb) = p(0 (pW , U W), w (p w , ) 


( 22 ) 


with pW and u W the density and the internal energy per unit mass 

f/W 


u (i) = 


m 


(0 


(23) 


The second term depends on the number of elements (n e ), the reference volume ( Vo ^) for 
element j, and the strain energy per unit volume in shear (V^), here assumed to be 


= (1 - d U) ) /i (i) tr (E e ^ T E e ^ 


(24) 


where d (j> is a shear damage variable and is a shear modulus. The third term depends 
on the number of subelements per element (n s ), the subelement reference volumes (Vo’®), 
and the strain energy per unit volume in tension (' 0 b») ; here assumed to be 


= - (1 - D®) K& < - 1 > 2 


(25) 


where D <3> is a normal damage variable, K (j> is a bulk modulus, < x > denotes the bracket 
function 

< x > = x u{x) (26) 

and u denotes the unit step function. Since the subelement Jacobians and the shear strain 
tensor depend on the particle center of mass coordinates 


= j(j,k) ( c (0) ) E e b) = E e(j) (c w , E ?j(j) ) 


it follows that the system potential energy has the general functional form 


(27) 


V = V (p<‘> , U {t) , c<‘> , d (i) , D u) , E p0 > ) 


(28) 
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The system potential energy defines the generalized conservative forces 


(i) P W _ dV 
pW 2 <9pW ’ 



dV 
dc W ’ 


as well as the deviatoric stress 

S 0') 

and the strain energy release rates 


1 dV 

y<j) dE p (?) 





9TW ’ 




dd& 


( 29 ) 


( 30 ) 


( 31 ) 


due to damage evolution. Note that when the internal energy is introduced as a generalized 
coordinate, the associated generalized forces are constant. 

With the system Lagrangian now defined, the next four sections describe evolution equa- 
tions for the internal state variables. 


DENSITY EVOLUTION RELATIONS 


This section derives density evolution relations for the particles, by extending certain 
exact results for the deformation kinematics of a unit cell of spherical particles, arranged 
in a body centered cubic packing scheme. For uniform compression of such a unit cell, in 
isolation, the cell center density pW is related to the reference density p 0 ^\ the center of 
mass separation distances and particle radii h (j) for its eight neighbors by 


pf 8 h 


f2(3h^\ 3 


( 32 ) 


If the particles are not spherical but ellipsoidal then 




?l = i + - y 

p« « 2 - 


3 = 1 - 


(Hd) 


( 33 ) 


This expression for an isolated cell may be extended to a cell array of arbitrary size by 
describing the same kinematics in rate form. Taking the time derivative of the last equation 
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and multiplying by a step function which allows for contact with near neighbors only 


yields 
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(0 


Pt 


(i) 


3 y 


3 = 1 




( 34 ) 


where the summation is now over all n particles. The coefficient bbbh) must satisfy 


0 if i—j 


w 


= < 0 if i^j and > % 


. (i) 

1 if i^j and C (i,j)3 < J§y 

in order to correctly reflect the dependence of particle contact distance on the local density. 
Otherwise particles will contact remote neighbors through intervening matter. Hence 


W^ j) = (1 - 8ij) u 1 - C (i,j) 
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( 35 ) 


where denotes the Kronecker delta. Note that W^ 1 ’ 31 performs no interpolation. Intro- 
ducing the kinematic relation for Q l,3 \ developed in an earlier section, yields 

3 n 
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( 36 ) 


which is the constraint form of the density evolution relations. The coefficients of the particle 
translational velocities and angular velocities in this expression will determine generalized 
forces in the momentum balance (Lagrange) equations derived in a later section. When the 
particle velocities in this expression are eliminated in favor of the particle momenta 


c « = m (i)- 1 p (i) , 




( 37 ) 


the density evolution equations take an explicit state space form convenient for use in nu- 
merical simulation. 


PLASTICITY AND DAMAGE MODELS 

This section introduces evolution equations for the plastic and damage variables. As in the 
case of the potential energy, alternative constitutive assumptions may be introduced without 
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change to the basic modeling methodology. The plastic flow rule used here is adapted from 
reference [27], and represents the simplest possible accommodation of the aforementioned 
isochoric plastic deformation constraint. The flow rule is 


]t )pU) — 


A(i) 


N p b) ]\f g p(j) 


||5 p b')|| = = 

where AbO is a positive proportionality coefficient, S p b* is an effective stress 

§ pU) — ]\j T f^pU) T g O') 


(38) 


(39) 


and the invariant operator is defined by 


ITII = 


* tr (T T T) 


1/2 


(40) 


for any second order tensor T. The fourth order tensor coefficients in the flow rule are 
defined by 


]\jpb) y _ 


1 


2 1 1 Ob') 1 1 
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]^T = T- -tr{ T)I 

for any symmetric second order tensor T. The yield function is 


(41) 

(42) 


f(j) = || j 5 'pC?)|| _ yb) 


(43) 


where Y (j > is the yield stress 


yb) = yb) (i _ d {j) ) (1 + K ^ e p{j) ) aW (l - r) U) 6 HU) ) 


(44) 


with e p b) the effective plastic strain, k ( j> a strain hardening coefficient, a b) a strain hardening 
exponent, 77 b) a thermal softening coefficient, and 9 H b) the homologous temperature. The 
effective plastic strain is determined by integrating the rate relation 


gPb) — ||E p bO| 


(45) 
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while the incremental plastic strain for a time step At is computed using 


*><,-) <||s»«||-y«> 

AA (l-d0))2/ i fa) 


( 46 ) 


The damage evolution equations applied here are adapted from reference [28], and dissi- 
pate the strain energy stored in tension and shear over fi time steps, once an element meets 
any stipulated material failure criteria. The evolution equations are 


D (*> = u( 1 - D U) ), 
hAt K h 


d U) = u( 1 - d U) ) 


(47) 


where is initialized to zero, and is set to a value of one when the accumulated plastic 
strain, temperature, or element compression reach corresponding critical values for the plastic 
failure strain (e^ ), melt temperature (dm), or maximum compression (JrP)- Other failure 
criteria may of course be specified. 

In general terms, the plastic and damage evolution equations are noholonomic constraints 
of the form 


E?(j) 

= E p(j) (p w , c w , d u \ D ij \ e p{j \ E p(j) ) 

(48) 

d U) 
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on the system level Lagrangian model. 


ARTIFICIAL VISCOSITY AND HEAT DIFFUSION 

Shock physics codes of the continuum or particle type incorporate a numerical viscosity 
and artificial heat diffusion. The forms used here are typical of particle codes, with one 
exception. Since the ellipsoidal particles used here admit rotational degrees of freedom, a 
viscous torque has been added which damps the relative rotation of neighboring particles. 

A viscous force is introduced for converging particles only 
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where the relative normal velocity is 
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and the viscosity coefficient is 
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with Cs'^ and V} 1 ' 1 a sonndspeed and particle reference volume. The parameters c Q and c\ are 
nondimensional linear and quadratic numerical viscosity coefficients. 

Similarly the viscous torque is 


M (i) = R (i)r ( R w w (i) - R^w^ ) u( 1 - 

3 = 1 


where the torsional damping coefficient is 
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Finally the thermal power flow due to artificial heat diffusion is taken to be 

n 

Q con{i) = R{i,j) ( 0{i) - dU) ) “(! - C (i,j) ) 

3 = 1 

where the heat transfer coefficient is 


(56) 


R W) = k f cfc^vy + 


(57) 


with ct a specific heat and k 0 a numerical heat diffusion coefficient. 


INTERNAL ENERGY EVOLUTION EQUATIONS 

The last internal state variable to be considered is the internal energy. The introduction 
of internal energy states as generalized coordinates allows the thermomechanical problem of 
interest here to be solved using energy methods. 

The internal energy evolution equations for particle i are 

lj(i) — jjwrk(i ) _|_ jjirr(i ) _ jjcon{i ) Pggh 
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where the first term accounts for mechanical work, the second term accounts for the effects 
of irreversible entropy production, and the third term represents numerical heat diffusion. 
The mechanical power flow for particle i is 


p{() 

t rwrk(i) _ (i) -(i) 

L ~ m ~(i) 2 P 


P 


(59) 


The energy dissipation due to irreversible entropy production for particle i depends on the 
viscous forces and torques, which act on the particles, and on the dissipation in the elements 

n 

jjirr(i) = f (i) T c W + M« T u> w + <j>M Q irr (60) 

3 = 1 

where Q irr(j > is a power flow due to damage evolution and plastic deformation in element j 


Q irr (J) = T d ^D U) + T dU) d U) + V Q e{j) tr (s {j)T E pU) ^j (61) 

and is the fraction of the dissipation in element j associated with particle i. 

Finally the internal energy flows due to numerical heat diffusion are 


^ jcon(i ) 
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(62) 


As in the case of the density evolution equations, a constraint form of the internal evolution 
relations is used to identify the generalized forces which appear in the Lagrange equations 
developed in the next section. For numerical implementation of the method, the generalized 
velocities are eliminated by introducing the momentum states as well as evolution relations 
for the density, plastic, and damage state variables. The resulting internal energy evolution 
relations take an explicit state space form. 


LAGRANGE’S EQUATIONS 

The preceding sections defined stored energy functions and nonholonomic constraints for 
the thermomechanical particle-element system. This section develops the final ODE model. 
The results of Shivarama and Fahrenthold [21] allow in the present case Lagrange’s equations 


to take the canonical form 
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where q c W, Q C W, Q p W, Q u< ^\ Q d ^\ Q D< ^\ and Q p V) are generalized forces determined by 
the nonholonomic constraints. The degenerate forms of the Lagrange equations for the 
internal state variables are due to the fact that those variables are not associated with any 
generalized momenta. Introducing Lagrange multipliers 7 p h) ; ^u(i) ^ y(j) , ^D{j) ^ arK ] 
for the constraints, the generalized forces are found to be 
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These results allow the unknown Lagrange multipliers to be determined in closed form, so 
that the final Lagrange equations are 
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where the generalized forces and torques due to particle interactions are 
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Supplemented by the evolution equations for density, internal energy, shear damage, normal 
damage, and plastic strain, the result is an explicit first order ODE model for the thermo- 
mechanical particle-element system. 


EXAMPLE SIMULATIONS 

This section describes four example problems which illustrate application of the improved 
particle-element method developed in this paper. The first two examples involve one dimen- 
sional test problems with known exact solutions. The third and fourth example problems 
involve three dimensional simulations. The third example models a published experiment 
while the fourth example measures the relative computational cost of the present method in 
a hypervelocity impact application of current research interest. 

The first example is the wall shock problem of Noh [31]. The simulations employ an ideal 
gas equation of state 

0 = — (u - u 0 ) 
c v 


P = (7 - 1) P (u - O, 


(80) 


and the parameters shown in Table 1. This problem models the collision of a fluid stream 
located initially in the region 0.0 < x < 0.5 with a rigid wall located at x — 0. The initial 
conditions are p = p D , u = u 0 , and v = —1. Figures 1 through 4 plot the simulation results for 
velocity, density, pressure, and temperature at the stop time of 0.3, for a model composed of 
200 particles. The numerical results show good agreement with the exact solution, although 
better results have been obtained using finite difference and finite element methods [31,35]. 
Table 2 shows convergence of the simulation results, as the particle count is increased, in 
terms of the velocity error norm 
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and the temperature error norm 


1/2 


6 0 = ' n 


\e w - 0(x w ) 


2=1 


(82) 


where v and 9 denote the exact solutions for the velocity and temperature. 

The second example is the bar impact problem of Kolsky [33]. The simulations employ a 
linear isothermal equation of state 


P=K (/>/ A> - 1) (83) 

and the parameters shown in Table 3. This problem models the one dimensional motion 
of an elastic bar of length L subjected to a step tensile pressure loading of magnitude P ex t, 
applied at the end x = L at time t = 0. Figures 5 and 6 plot the simulation results for the 
bar midpoint velocity versus time, to a stop time of 0.001 seconds, for two different particle 
counts. The numerical results show good agreement with the exact solution. The results 
shown in Figure 5 are quite similar to those reported by Lu et a. [34], at the same particle 
count, employing an element-free Galerkin method with an explicit integration algorithm. 

The third example problem models the oblique impact of a tungsten alloy (DX2HCMF) 
rod on a steel (SIS 2541) plate, an experiment described in references [29] and [30]. The 


simulations employed a Mie-Gruneisen equation of state [32] and the material properties 
[30,32] listed in Table 4. The cylindrical projectile has a diameter of 0.5 cm and a length 
of 7.5 cm (L/D = 15). The simulation models a 1.5 krn/s impact on a 0.5 cm thick plate 
at a sixty degree obliquity, and was run at three different particle counts. Figures 7 and 8 
show the initial configuration and the simulation results at 100 microseconds after impact, 
while Figures 9 and 10 show sectioned views at 20 and 40 microseconds after impact. This 
simulation illustrates the fact that the present method retains all material fragments and 
models contact-impact of all intact and fragmented material. Table 5 provides simulation 
results, at several particle counts, for the residual rod length and residual rod velocity, 
showing good agreement with the corresponding experimental values (6.38 cm and 1.46 
km/s). The simulation results shown here differ by less than three percent from those 
reported by Lee and Yoo [30] for Lagrangian finite element methods. 

The fourth example problem models the oblique impact of an aluminum sphere on a re- 
inforced carbon-carbon plate. This example does not model a specific experiment; rather it 
was used to measure the improved computational efficiency of the method described here, 
in a hypervelocity impact application of current research interest [6,36]. The projectile ma- 
terial, diameter (0.618 cm), and impact velocity (7 km/s) represent a typical orbital debris 
impact threat [4], while the target plate thickness (0.47 cm) is characteristic of reinforced 
carbon-carbon components of the Space Shuttle thermal protection system [6]. Experimen- 
tal studies of the mechanical properties of reinforced carbon-carbon are in progress; the 
material properties used in the simulations are listed in Table 6. A total of six simulations 
were performed, at three different particle counts, to measure the the computational cost 
of density calculations made using interpolations kernels, as compared to the nonholonomic 
formulation developed in the present paper. Each simulation was run to a stop time of ten 
microseconds, sufficient to perforate the target and include the shock loading and fragmen- 
tation processes of central interest in hypervelocity impact applications. Figures 10 through 
15 show representative simulation results. Table 7 lists the wall clock times and processor 
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counts for the simulations, run on a Linux cluster composed of dual processor (3 GHz) nodes. 
The results show the relatively high computational cost of kernel based density calculations, 
in the present hybrid particle finite element context. Wall clock times are increased for all 
three problem sizes and associated processor counts, by an average factor of one third. This 
result is not surprising, since the finite element related portion of the computation is rela- 
tively inexpensive [22] , while most of the particle related computations are performed in two 
sequential routines: one loops over all neighbor particles to determine the density, the second 
loops again over all neighbor particles to compute the particle interaction forces. The effect 
of introducing the nonholonomic density calculation developed here is to eliminate the first 
of these two routines. Considering the very high computational costs of three dimensional 
shock physics problems, the measured reduction in wall clock time is significant. 

CONCLUSION 

The present paper has formulated a new kernel free particle-finite element method and 
demonstrated its application in three dimensional hypervelocity impact simulations. Unlike 
alternative methods, the formulation is derived without reference to any weighting or in- 
terpolation functions for either the density or the rate of dilatation. The improved method 
introduces a new formulation of the thermomechanical Lagrange equations, one which em- 
ploys internal energies as generalized coordinates. This avoids the requirement to perform 
certain Legendre transforms, in order to express the dependence of the pressure and tem- 
perature on entropy, and hence allows for the use of equations of state in standard form. As 
compared to the previous formulation, the revised method is both simplified and computa- 
tionally more efficient. Applications work in progress is focused on the simulation of orbital 
debris impact effects on spacecraft thermal protection materials [6]. 
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Table 1. Simulation parameters for the wall shock problem 


Ratio of specific heats ( 7 ) 

5 

3 

Specific heat (c„) 

1.0 

Reference internal energy (e D ) 

1.0 

Reference density (p Q ) 

1.0 

Reference temperature (0 o ) 

1.0 

Numerical viscosity coefficient (c 0 ) 

1.0 

Numerical viscosity coefficient (ci) 

0.0 

Numerical conduction coefficient (k 0 ) 

0.1 


Table 2. Error norms for the wall shock problem 


Number of particles 

Velocity error norm 

Temperature error norm 

100 

0.17116 

0.08594 

200 

0.12141 

0.06075 

400 

0.08635 

0.04325 

800 

0.06081 

0.03032 


Tabic 3. Simulation parameters for the bar impact problem 


Length of the bar (A, in) 

100 

Bulk modulus (A", psi) 

30.0 x 10 6 

Applied end loading ( P ext , psi) 

50.0 x 10 3 

Reference density (p D , pci) 

0.73 x 10 “ 3 

Numerical viscosity coefficient (c G ) 

1.0 

Numerical viscosity coefficient (ci) 

0.0 
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Table 4. Material properties used in the long rod impact simulations 


Material property 

Projectile 

Target 

Reference density (g/cc) 

17.6 

7.87 

Shear modulus (Mbar) 

1.45 

0.801 

Reference yield stress (Mbar) 

0.0075 

0.0105 

Strain hardening coefficient 

1.15 

0.177 

Strain hardening exponent 

0.49 

0.12 

Thermal softening coefficient 

1.0 

1.0 

Melt temperature (deg K) 

1,700 

1,723 

Specific heat (Mbar-cm 3 per g-deg K) 

0.143e-5 

0.448e-5 

Plastic failure strain 

1.0 

1.0 


Table 5. Simulation results for the long rod impact problem 


Number of particles 

Residual length 
(cm) 

Residual velocity 
(km/s) 

35,992 

6.08 

1.42 

92,498 

6.47 

1.42 

187,826 

6.59 

1.42 


Table 6. Material properties used in the plate impact simulations 


Material property 

Projectile 

Target 

Reference density (g/cc) 

2.70 

1.58 

Shear modulus (Mbar) 

0.271 

0.0718 

Reference yield stress (Mbar) 

0.0029 

0.000771 

Strain hardening coefficient 

125.0 

2.0 

Strain hardening exponent 

0.10 

1.0 

Thermal softening coefficient 

0.567 

-1.0 

Melt temperature (deg K) 

1,220 

3,840 

Specific heat (Mbar-cm 3 per g-deg K) 

0.884e-5 

0.712e-5 

Plastic failure strain 

1.0 

0.5 
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Table 7. Relative computational costs for the plate impact simulations 


Particles 

Density calculation 

Processors 

Wall clock hours 

Relative cost 

21,334 

nonholonomic 

4 

1.2611 

1.000 

21,334 

kernel 

4 

1.7581 

1.394 

63,253 

nonholonomic 

6 

4.5047 

1.000 

63,253 

kernel 

6 

5.8078 

1.289 

140,070 

nonholonomic 

8 

28.8933 

1.000 

140,070 

kernel 

8 

38.1411 

1.320 
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Figure 1. Wall shock problem, velocity versus position at t = 0.3. 

Figure 2. Wall shock problem, density versus position at t = 0.3. 

Figure 3. Wall shock problem, pressure versus position at t = 0.3. 

Figure 4. Wall shock problem, temperature versus position at t = 0.3. 

Figure 5. Bar impact problem, midpoint velocity versus time for 51 particles. 

Figure 6. Bar impact problem, midpoint velocity versus time for 101 particles. 

Figure 7. Long rod impact problem, element plot of the initial configuration 

Figure 8. Long rod impact problem, particle-element plot at 100 microseconds after 
impact 

Figure 9. Long rod impact problem, sectioned particle-element plot at 20 microseconds 
after impact, color on temperature. 

Figure 10. Long rod impact problem, sectioned particle-element plot at 20 microsec- 
onds after impact, color on temperature. 

Figure 11. Plate impact problem, element plot of the initial configuration. 

Figure 12. Plate impact problem, particle-element plot at 10 microseconds after im- 
pact. 

Figure 13. Plate impact problem, element plot at 10 microseconds after impact. 

Figure 14. Plate impact problem, sectioned particle-element plot at 10 microseconds 
after impact, color on temperature. 

Figure 15. Plate impact problem, sectioned element plot at 10 microseconds after 
impact, color on effective plastic strain. 
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Figure 1: Wall shock problem, velocity versus position at t = 0.3. 
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Figure 2: Wall shock problem, density versus position at t = 0.3. 
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Figure 3: Wall shock problem, pressure versus position at t = 0.3. 
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Figure 4: Wall shock problem, temperature versus position at t = 0.3. 
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midpoint velocity versus time 



Figure 5: Bar impact problem, midpoint velocity versus time for 51 particles. 
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Figure 6: Bar impact problem, midpoint velocity versus time for 101 particles. 
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Figure 7: Long rod impact problem, element plot of the initial configuration. 
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Figure 9: Long rod impact problem, sectioned particle-element plot at 20 microseconds after 
impact, color on temperature. 
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Figure 10: Long rod impact problem, sectioned particle-element plot at 40 microseconds 
after impact, color on temperature. 
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Figure 11: Plate impact problem, element plot of the initial configuration. 
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Figure 12: Plate impact problem, particle-element plot at 10 microseconds after impact. 
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Figure 13: Plate impact problem, element plot at 10 microseconds after impact. 
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Figure 14: Plate impact problem, sectioned particle-element plot at 10 microseconds after 
impact, color on temperature. 
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Figure 15: Plate impact problem, sectioned element plot at 10 microseconds after impact, 
color on effective plastic strain. 
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INTRODUCTION 

Advanced orbital debris shield designs often incorporate a multi-plate geometry and are 
fabricated of composite materials. An example is the aluminum-Nextel-Kevlar shields de- 
veloped by NASA [1] for application on the International Space Station (ISS). Simulation 
of hypervelocity impact effects on these shielding designs is difficult, due both to their ge- 
ometry and their material composition. Since multi-plate shields distribute the allocated 
shielding mass over more than one bumper, they reduce the minimum bumper thickness and 
in general increase the computational cost of the simulation. The extensive use of composite 
materials in advanced shielding designs complicates the simulation problem for two basic 
reasons: (a) the requirement to consider nonhomogeneous or anisotropic media increases the 
number of internal state variables which must be evolved in any numerical simulation, and 
(b) the thermomechanical response of such materials, generally more complex than that of 
aluminum or other metals, may be only partially described by the existing material property 
data base. In spite of these complications, simulation of the hypervelocity impact response 
of multi-plate metal-composite debris shields is of strong engineering interest, due to their 
proven effectiveness as components of spacecraft protection systems [2], The present work 
describes numerical implementation of a rate dependent material model for Kevlar [3,4], and 
the conduct of several three dimensional simulations, performed to investigate the hyper- 

3 Professor, corresponding author, phone: (512) 471-3064, email: epfahren@mail.utexas.edu 

4 Graduate research assistant 
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velocity impact response of aluminum-Nextel-Kevlar debris shields like those deployed on 
the ISS. The simulations employed a hybrid-particle finite element method and a parallel 
computer code specifically developed to address the orbital debris shielding design problem. 
Simulation work appears to conservatively estimate the protection afforded by multi-layered 
aluminum-Nextel-Kevlar shielding. 

COMPOSITE DEBRIS SHIELDING 

Kevlar aramid fiber, introduced by the Du Pont Company in 1972, has been used in a 
variety of impact protection applications. The woven cloth Nextel (manufactured by 3M 
corporation) is made from alumina, a ceramic shown to provide effective impact protection 
when applied as a component of composite armor systems. These materials are used in 
the design of the multiplate orbital debris shielding deployed on some modules of the ISS. 
The latter shielding consists of three material layers, arranged in the sequence: aluminum, 
Nextel, Kevlar. The outer aluminum bumper is located approximately eleven centimeters 
from the pressurized module, a multilayer stack of Nextel cloth is located approximately 
halfway between the aluminum bumper and the pressure wall, and a multilayer stack of 
Kevlar cloth is located directly behind the Nextel. Nextel is manufactured in various grades. 
A number of different grades of Kevlar are also manufactured, with Kevlar 29 (all purpose 
yarn), Kevlar 49 (high modulus yarn), and Kevlar 129 (high tenacity yarn) used in composite 
armor applications. Tables 1 and 2 list mechanical properties of these three grades of Kevlar 
and compare some important material properties of high performance Nextel and Kevlar 
fibers. 

IMPACT SIMULATIONS 

Experimental work has demonstrated that aluminum-Nextel-Kevlar shields perform much 
better than weight equivalent Whipple shield designs. Hence the development of a validated 
computer aided design tool for use in future composite shielding design studies is of consid- 
erable interest. Fahrenthold and Shivarama [9] described some initial work on the use of a 
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hybrid particle-finite element method for composite shield impact simulations. The latter 
work used a rather coarse three dimensional model of less the 0.5 million particles, and a 
simulation time of less than 55 microseconds. The simulations discussed here involved an 
order of magnitude increase in particle count and a three-fold increase in simulation time. 
Related work using SPH and finite element methods has been reported by Hiermaier et al. 
[10] and Palmieri [11], although the latter efforts have focused on the use of Kevlar-epoxy 
plates [12], in lieu of the cloth Kevlar material modeled here. Three of the four simula- 
tions described in the present chapter modeled an experiment performed by Groscli [13], 
test number SwRI 7139-24, in which a 1.07 gram inhibited shaped charge (ISC) projectile 
struck a two-thirds scale model composite orbital debris shield similar to that deployed on 
the ISS. Parameters of the experiment are detailed in Table 3. The impact velocity is 11 
kilometers per second, and the impact obliquity is 45 degrees. The first simulation reported 
here modeled the ISC projectile used in the experiment. A second simulation was performed 
using a spherical projectile with the same mass as the experimental (ISC) projectile, in or- 
der to investigate the effect of projectile shape [14] on the simulation results. A concern in 
the conduct of any hypervelocity impact simulation is the use of an appropriate equation 
of state. However there is no validated Mie-Gruneisen or tabular equation of state data 
for either Nextel cloth or Kevlar cloth. The aforementioned work of Hiermaier et al. [10] 
developed an equation of state and an orthotropic elastic model for Kevlar-epoxy, as well as 
a compaction equation of state for Nextel. The Kevlar layer modeled here does not include 
epoxy, and a compaction model for Nextel may not be appropriate for use in the present 
application. Hence the first and second simulations presented here used a linear equation of 
state and the material properties listed in Table 4. The first simulation represented the ISC 
projectile as a hollow aluminum cylinder, with dimensions estimated from the experimen- 
tal radiographs. The exact mass and geometry of ISC projectiles is somewhat uncertain. 
Both simulations employed approximately 6 million particles, and required approximately 
252 wall clock hours to complete, in parallel execution on 16 processors of an IBM Regatta. 
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The enclosed Figures 1 and 2 show an element plot of the first case at the simulation start 
time, and a particle plot of the simulation results at 100 microseconds after impact. Figures 
3 and 4 show element plots of the damage at each level of the multi-plate structure, again at 
100 microseconds. The simulation results underestimate the composite shield performance, 
since the experimental results showed a bulged but not perforated wall plate. A comparison 
of the first and second simulations, involving mass equivalent spherical and ISC projectiles, 
is provided by the wall plate plot of Figure 5, showing the results of the spherical projectile 
impact simulation, at 150 microseconds after impact. The results suggest that projectile 
geometry effects are significant even at the upper end of the orbital debris velocity regime. 

COMPOSITE MATERIAL MODELS 

The preceding simulation results underestimate the protection afforded by composite 
shielding materials. To investigate the effects of changes in the composite material mod- 
els on the predicted shielding performance, both the Nextel and Kevlar equations of state 
were modified, and a rate-dependent strength model was introduced for the Kevlar. Two 
additional simulations were then conducted using the modified composite material models. 
The modified equations of state for Nextel and Kevlar used a Mie-Gruneisen functional form 
and the estimated material properties listed in the enclosed table, with the Gruneisen gamma 
for Nextel and Kevlar approximated using the formula of Anderson [17] 

T = (3k/ P C (1) 

where (3 is the thermal expansion coefficient, k is the bulk modulus, C is the specific heat, 
and p is the density. The engineering application discussed here involves very high strain 
rate loading, hence any variation of material properties with strain rate is of considerable 
interest. Two recent papers [3,4] describe a significant strain rate dependence in the measured 
mechanical response of Kevlar 49, as summarized in Table 5, for experiments conducted over 
a large strain rate regime. The parameters listed refer to the (apparent) Young’s modulus, 
the maximum engineering stress, and the engineering strain at maximum engineering stress, 
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as measured under uniaxial loading conditions. There is no data showing a similar strain 
rate dependence for Nextel. For use in later impact simulation work, the maximum stress 
versus strain rate data were fit to the function 

- = 9(e) (2) 

0o 

where 


g(e) = 1 i < £ 0 


and 

9(e) = 1 + a 

As indicated in Figure 6, for Kevlar 49 


In 

eo 


e > £q 


a = 0.0064, m = 1.4 
where the reference values of stress and strain rate are 


( 3 ) 


( 4 ) 


( 5 ) 


a 0 = 2340 MPa, i 0 = 0.0001 s” 1 (6) 

To investigate the effects of rate dependence on predicted Kevlar impact response, the 
plasticity model of Fahrenthold and Horbon [5] was modified, as described in the equations 
which follow, and used in the two simulations discussed next. Note that the plasticity 
model outlined here incorporates large strain kinematics and an isochoric plastic deformation 
constraint. The nonassociated flow rule for the plastic strain rate (E P ) is 

E p = A (C P W + WC P ) (7) 

where A is a scalar multiplier and 

w = C p S e// + S e// C p - —tr [C p S e// + S e// C p ] I (8) 

3 
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with the effective stress tensor and plastic Cauchy-Green strain tensor defined by 


S e// = - — d)2fj, 0 E e ^ C p = 1 + 2E p ( 9 ) 

9\£) 

where e is the effective strain rate, E e is the elastic deviatoric strain, d is the deviatoric 
damage, and /i 0 is the shear modulus. The yield condition is 


f = r-Y, 


T = 



ge// T ge// 



( 10 ) 


where r is the second invariant of the effective stress. The yield stress Y is 


r = i(i-d)io {i-7»} 


(ii) 


where Yq is the reference yield stress, 7 is a thermal softening modulus, and 6 is the homolo- 
gous temperature. The plastic strain increment at each time step is determined using a one 
step iteration procedure with 


M (r - Y) A (r - Y) 
(1 - d) 2/i 0 u :g(i) 

where A denotes the unit step function and 

ui = jtfr [W T W] } 


( 12 ) 


(13) 


The preceding modified material models were used in two additional three dimensional 
simulations of hypervelocity impact effects on aluminum-Nextel- Kevlar shielding. The first 
simulated a published light gas gun experiment [1], while the second again considered the 
ISC experiment discussed in the last section. The light gas gun experiment was JSC test 
B536, and involved the oblique (15 degree) impact of a 1.0 gram aluminum sphere on a full 
scale aluminum-Nextel- Kevlar shield, at a velocity of 6.86 krn/s. Parameters of the experi- 
ment are listed in Table 5. The simulation employed 1.18 million particles and required 63.2 
wall clock hours to complete on 16 processors of an IBM Regatta. The simulation results 
(Figures 7 and 8) indicate at 150 microseconds a slightly deformed but not perforated wall 
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plate. The experiment showed a “slight dish,” indicating that the simulation results provide 
an accurate estimate of the wall plate damage. Since simulation of this light gas gun test 
showed good agreement with the corresponding experiment, simulation of the ISC projectile 
test SwRI 7139-24 was repeated with a relatively high resolution model and the modified 
material models discussed in this section. The simulation employed 15.6 million particles 
and required 305 wall clock hours to complete on 512 processors of an SGI Origin. The sim- 
ulation results (Figures 9 and 10) indicate at 75 microseconds a perforated wall plate. The 
effect of introducing nonlinear equations of state for the composites and a rate dependent 
Kevlar model appears in this case to be rather small, with the simulation again overestimat- 
ing the wall plate damage. The results reported here suggest that the hybrid particle-element 
method used in the simulations is numerically robust and captures important basic features 
of multi-plate impact experiments. For example, the intermediate composite shields are 
perforated, not fluidized, by the debris cloud impact. However the results also suggest that 
modeling improvements are needed, in order to better represent the available experimental 
impact data. The present study suggests possibilities for future work in both the model- 
ing and experimental areas. First, the modeling of damage induced anisotropy, commonly 
encountered for example in ceramic media, may be investigated to determine its effect on 
predicted shield performance. Second, impact experiments using a simpler target geometry 
might be conducted, with a primary goal of materials characterization in mind. 

CONCLUSION 

The hybrid numerical method utilized here has been applied with some success to model 
the impact of metal projectiles and targets, over a wide range of velocities. The present 
work has described development of a rate dependent material model for Kevlar and three 
dimensional simulation work aimed at extending the formulation, for use in the computer 
aided design of multi-plate composite orbital debris shields. Both the geometry and the ma- 
terial composition of such shields complicate the simulation problem. The simulation results 
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conservatively estimate the protection afforded by multi-layered aluminum-Nextel-Kevlar 
shielding. Additional material modeling work, guided by impact tests focused on materials 
characterization issues, is expected to improve upon current simulation capabilities. The 
rather high computational cost of three dimensional multi-plate shield impact simulations 
emphasizes the importance of developing efficient parallel numerical implementations. 
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Yarn properties 

Kevlar 29 

Kevlar 49 

Kevlar 129 

Tensile Strength (Kpsi) 

420 

420 

485 

Initial Modulus 

10.3 

17.4 

14 

Elongation (%) 

3.6 

2.8 

3.3 

Density (lb/ft 3 ) 

89.9 

90.5 

90.5 


Table 1: Mechanical properties of Kevlar ararnid fibers 


Fiber type 

Density 

(g/cm 3 ) 

Strength 

(GPa) 

Modulus 

(GPa) 

Elongation 

(%) 

Fiber diameter 

w 

Kevlar 129 

1.45 

3.4 

99 

3.3 

12 

Nextel 

2.5 

1.72 

152 

2 

13 


Table 2: Comparison of mechanical properties for Kevlar 129 and Nextel 


Projectile mass (aluminum, L/D = 1) 

1.07 g 

Bumper thickness (aluminum) 

0.127 cm 

Nextel areal density 

0.400 g/cm 2 

Kevlar areal density 

0.128 g/cm 2 

Wall plate thickness (aluminum) 

0.3175 cm 

Total standoff 

7.62 cm 

Projectile velocity 

11.25 km/s 

Impact obliquity 

45 degrees 


Table 3: Parameters of SwRI test number 7139-24 
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Material property 

Aluminum 

Nextel 

Kevlar 

Shear modulus (Mbar) 

0.271 

0.164 

0.100 

Reference density (g/cc) 

2.7 

2.7 

1.45 

Reference sound speed (cm//isec) 

0.524 

0.4968 

0.5352 

Reference yield stress (Mbar) 

0.0029 

0.0172 

0.034 

Strain hardening exponent 

0.1 

0.0 

0.0 

Strain hardening modulus 

125.0 

0.0 

0.0 

Melt temperature (degrees Kelvin) 

1,220 

1,220 

700 

Specific heat 




(Mbar-cm 3 per g-kilodegrees Kelvin) 

0.00884 

0.00884 

0.00142 

Spall stress (Mbar) 

0.012 

0.100 

0.100 

Plastic failure strain 

1.0 

1.0 

1.0 

Thermal expansion coefficient 




(per kilodegrees Kelvin) 

0.0216 

0.009 

0.038 

Mie-Gruneisen gamma 

1.97 

0.2513 

0.7666 

Mie-Gruneisen slope coefficient 

1.4 

1.0 

1.0 


Tabic 4: Material properties used in the simulation 


Strain rate 

0.0001 s- 1 

0.01 s" 1 

140 s — 1 

440 s — 1 

1350 s” 1 

E (GPa) 

97 

100 

112 

119 

125 

®max (GPa) 

2.34 

2.47 

2.94 

3.02 

3.08 

£m (%) 

3.29 

3.33 

3.54 

3.64 

3.86 


Table 5: Mechanical properties of Kevlar 49 versus strain rate 


Projectile mass (aluminum sphere) 

1.0 g 

Bumper thickness (aluminum) 

0.16 cm 

Nextel areal density 

0.600 g/ern 2 

Kevlar areal density 

0.192 g/cm 2 

Wall plate thickness (aluminum) 

0.48 cm 

Total standoff 

11.4 cm 

Projectile velocity 

6.86 krn/s 

Impact obliquity 

15 degrees 


Table 6: Parameters of JSC test number B536 
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Figure 1: Element plot of the initial configuration 
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Figure 2: 


Particle plot of the simulation results at 100 microseconds after impact 
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Figure 3: Element plot of the simulation results at 100 microseconds after impact 
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Figure 4: Element plot of the Kevlar shield and wall plate at 100 microseconds after impact 
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Figure 5: Element plot of the wall plate damage, spherical projectile 
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Figure 6: Strain rate dependence of tensile strength of Kevlar 49 
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Figure 7: Element plot of the initial configuration 
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Figure 8: Element plot of the Kevlar shield and wall plate at 150 microseconds after impact 
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Figure 9: Element plot of the simulation results at 75 microseconds after impact (front view) 
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Figure 10: Element plot of the simulation results at 75 microseconds after impact (rear view) 
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ABSTRACT 


Spacecraft operating in low earth orbit face a significant orbital debris impact hazard. 
Of particular concern, in the case of the Space Shuttle, are impacts on critical components 
of the thermal protection system. Recent research has formulated a new material model 
of reinforced carbon-carbon, for use in the analysis of hypervelocity impact effects on the 
Space Shuttle wing leading edge. The material model has been validated in simulations 
of published impact experiments and applied to model orbital debris impacts at velocities 
beyond the range of current experimental methods. The results suggest that momentum 
scaling may be used to extrapolate the available experimental data base, in order to predict 
the size of wing leading edge perforations at impact velocities as high as 13 krn/s. 

NOMENCLATURE 


d, shear damage variable 
D , projectile diameter 
D p , perforation diameter 
D c: coating spall diameter 

D. rate of deformation tensor 

e, Euler parameter vector 

E. deviatoric strain tensor 
E e , elastic strain tensor 
E p , plastic strain tensor 
e p , effective plastic strain 
e, deviatoric strain rate 

F. deformation gradient tensor 


Sij, Kronecker delta 
L. velocity gradient tensor 
N e , number of elements 

R, rotation matrix 

S, deviatoric stress tensor 
S p , effective stress tensor 
s, entropy density 

v. impact velocity 
Y, yield stress 
•0, strain energy density 
9, temperature 
< i i>, impact obliquity 
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INTRODUCTION 


Carbon-carbon composites offer an unusual combination of thermal and mechanical prop- 
erties. 1 Their light weight and high temperature strength satisfy some very stringent design 
requirements for reusable orbital vehicles. 2 The wing leading edge of the Space Shuttle, 
subject to severe thermal re-entry loads, is constructed of reinforced carbon-carbon (RCC) 
panels, coated in silicon carbide to prevent oxidation. 3 Although the thermal properties of 
RCC composites are well understood, 4 much less is known about their dynamic mechanical 
properties. The loss of the Space Shuttle Columbia, 5 apparently due to impact damage on 
the wing leading edge, has motivated recent experimental 6 and computational' work aimed 
at developing a better understanding of the impact response of thermal protection materials. 

The wing leading edge damage to Columbia was unexpected, the result of a relatively 
low velocity impact by a relatively low density projectile. 8 Another impact damage hazard, 
due to space debris in low earth orbit, has long been recognized. This threat involves 
projectiles of very low mass, but much higher density, and impact velocities as high 15 km/s. 
The debris shielding on the International Space Station is designed to defeat centimeter 
sized aluminum projectiles. Although the likelihood of such a projectile striking the Space 
Shuttle is quite low, orbital debris damage by much smaller projectiles is routinely observed 
during post-mission inspections of the vehicle. As a result previous experimental research 
has investigated the response of Space Shuttle thermal protection materials to orbital debris 
impact by spherical aluminum projectiles as large as 0.628 cm in diameter. 9 

Due to the high cost of carbon-carbon composites and the long fabrication lead times 
associated with the preparation of test samples, impact testing of RCC materials has been 
limited. In addition, the limitations of current experimental technology preclude hyperveloc- 
ity impact testing over the entire projectile mass and kinetic energy range of interest. As a 
result, numerical simulation can serve as an important complement to experimental studies 
of the impact response of RCC materials. Numerical models validated by comparison with 
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experiment at velocities below 8 krn/s can be used to extrapolate results into a higher veloc- 
ity impact regime. A coordinated experimental and computational approach to the study of 
RCC response to insulating foam impacts has proven to be productive; 10 the present paper 
extends the latter computational work, to projectiles and impact velocities associated with 
orbital debris impact. In particular it develops a new anisotropic, rate-dependent material 
model for reinforced carbon-carbon, validates that model in three dimensional simulations of 
published hypervelocity impact experiments, and applies the validated formulation in simu- 
lations of impacts at velocities beyond the experimental range. The results indicate that a 
momentum scaling approach used to correlate the available experimental impact data may 
be extrapolated to describe RCC perforation by hypervelocity projectiles at velocities as 
high as 13 krn/s. 

The present paper is organized as follows. The first section outlines the hybrid particle- 
finite element method used in the present study, including the imbedded large deformation 
kinematics and general functional forms for the associated constitutive relations. The second 
section discusses published experimental results on the properties of RCC. The third section 
develops an RCC constitutive model, formulated for use in hypervelocity impact applica- 
tions and reflecting important mechanical characteristics described in the material testing 
literature. The fourth section validates and applies the developed model in a series of three 
dimensional impact simulations. The last section presents conclusions and suggestions for 
related future work. 


NUMERICAL METHOD 


The material model described in this paper was developed for application in a specific 
numerical framework, the hybrid particle-finite element formulation of references 11 and 
12. In order to provide appropriate context, this section summarizes the latter numerical 
formulation, details certain element level kinematics, and provides functional forms for the 
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required constitutive relations. The kinematic and constitutive modeling framework assumed 
here has wide scope, so that the material model described in the present paper may be 
adapted for use in shock physics codes which are based on alternative numerical modeling 
schemes . 13 

The hybrid particle-finite element model employed here takes an explicit state space form. 
The state equations consist of evolution equations for the following variables: 

• translational and rotational momentum vectors for the three dimensional motion of 
ellipsoidal particles, 

• center of mass position vectors and Euler parameters for the particles, the latter pro- 
viding a singularity free description of particle rotations, 

• density and entropy for each particle, and 

• damage and plastic internal state variables for each finite element. 

The state equations are derived using a thermomechanical formulation of the Lagrange equa- 
tions. All inertia effects are modeled using the particles, whose mass centers are also nodal 
coordinates for the finite elements. The volumetric thermomechanical response of the mod- 
eled medium is described by an equation of state for the particles, which may take either an 
analytic or tabular form. 

The material modeling work described in the present paper develops two specific compo- 
nents of the general numerical formulation: 

• a strain energy density in shear, one part of the thermomechanical Lagrangian for the 
modeled particle-element system, and 

• a plasticity model which specifies evolution equations for the plastic internal state 
variables, equations which serve as nonholonomic constraints on the system level model. 
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The strain energy density in shear takes the general functional form 


F = if>(d,e,E,W) (1) 

where d is a shear damage variable, E is the total deviatoric strain, E p is the plastic strain, 
and e is a vector of Euler parameters which relates a material reference frame for each 
element to a single global Cartesian reference frame. The evolution equations for the plastic 
strain components take the general functional form 

E p = E p (s, d, e p , e, J, e, E. E p ) (2) 

where s is an entropy density, e p is the effective plastic strain, e is a deviatoric strain rate, 
and 

J = det{ F) (3) 

where F is the deformation gradient tensor. 

The strain and strain rate variables which appear in the preceding functional forms are 
defined by the following large deformation kinematics. 14 The deviatoric strain is 

E=t(C-l) (4) 

where 

C = F i F. F = (detF)-3F (5) 

The elastic shear strain is defined as 

E e = E - E p (6) 

where the flow rule for the plastic stain tensor must satisfy the isochoric plastic deformation 
constraint 

tr ( C P ~ T & ) =0, C p = I + 2 E p (7) 
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The effective plastic strain is determined by integrating the rate relation 


i p = 



( 8 ) 


with the indicated invariant operator defined by 

1 


IT 1 1 = 


- tr (T T T) 


1/2 


(9) 


for any second order tensor T. The deviatoric strain rate is 


e = | |D'| |, 

where D is the rate of deformation tensor 

1 


D ; = D — -tr( D) I 


( 10 ) 


D 




L = FF 


-l 


( 11 ) 


with L the velocity gradient tensor. 

In the case of anisotropic materials, the constitutive response is described in a material 
reference frame. Here an Euler parameter vector 


e = [e 0 ei e 2 e 3 ] r , e T e = 1 (12) 

is used to define a rotation matrix (R) for each element 



R 

= AG 

T 


(13) 


—ei 

e 0 - 

-e 3 

e 2 
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— e 2 
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G = 

-e 2 

— e 3 

e 0 

ei 

(15) 


. -e 3 

e 2 - 

-ei 

e o 



which relates a material reference frame in each element to the global Cartesian system used 
in the numerical simulations. The rotation matrix relates vector components p in the global 
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coordinate system to corresponding components q described in the material reference frame, 
using 

p = R q (16) 

The corresponding transformation relation for second order tensors is 

P = RQR t (17) 

where the components P refer to the global frame and the components Q refer to the material 
reference frame. 


REINFORCED CARBON-CARBON 


The published material property data base for carbon-carbon composites is limited by 
material costs and proprietary considerations. On the other hand, the complex nature of 
both the material and the application of interest here means that a rather wide range of 
experiments are needed to fully characterize its constitutive response. This section discusses 
some properties of RCC of particular significance in hypervelocity impact applications. 

The most directly relevant experimental results are those of Lu et ah, 6 who performed 
tests at Sandia National Laboratories on samples taken from Space Shuttle wing leading edge 
panels, in support of the Columbia accident investigation. They provide data on elastic mod- 
uli as well as strength measurements obtained from tension, bending, and compression tests. 
Although the elastic moduli measured in tension and compression were similar, strength in 
compression was approximately double that in tension. In addition they reported a strain 
rate dependence of the tensile strength, observing a fifteen percent increase in strength as 
the loading rate increased from 1 to 200 sec -1 . Finally they noted that removal of the silicon 
carbide coating from the tested samples showed little effect on the measured mechanical 
properties. 
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Several different authors have reported results of shear tests performed on carbon-carbon 
composites . 15-20 In the case of the RCC, interlaminar shear strength and stiffness is of 
interest, since oblique hypervelocity impacts will in general lead to multi-axial loading. The 
published shear test data show that interlaminar shear stiffness and strength can differ by 
factors of approximately two and four respectively, from their in-plane counterparts. 

Perhaps the most unusual property of RCC is its increase in strength with tempera- 
ture , 17,21 by as much as a factor of two, as compared to the thermal softening response 
observed in metals. The high temperature strength of RCC is important in hypervelocity 
impact applications, due to the adiabatic heating typically associated with shock loading. 

The preceding references, along with the equation of state literature 22 and published data 
on the thermal properties of RCC , 4,23 were used to estimate the material parameters used 
in the simulations reported in a later section. Although additional tests on Space Shuttle 
wing leading edge panels, like those reported by Lu at ah, are needed, the cited references 
represent the best data available to the authors at the time the simulations were conducted. 

MATERIAL MODEL 

Composite materials are used in structural , 24 orbital debris shielding , 25 and thermal 
protection 4 applications on a variety of spacecraft, hence their response to hypervelocity 
impact effects has been analyzed in a number of previous experimental 26 and computa- 
tional 27 studies. Previous material modeling work has considered both micromechanical 28 
and anisotropic continuum models. The present paper employs an anisotropic continuum 
approach, since the high computational cost of micromechanical models normally precludes 
their use in structural scale simulations. Large deformation, anisotropic continuum models of 
composite materials 29,30 normally address shock physics problems by extending small strain 
formulations originally developed for applications in structural mechanics. In an alternative 
approach, the present work starts with the finite strain, hybrid particle-element kinematics 


144 


discussed in an earlier section, and then formulates: (1) an anisotropic strain energy den- 
sity function which depends on a general deviatoric Lagrangian strain tensor, and (2) an 
anisotropic, temperature and rate dependent plastic flow rule which depends on an effective 
stress 31 and satisfies a general isochoric deformation constraint. Both the strain energy den- 
sity function and the plastic flow rule: (1) account for differences in material response under 
tension and compression, (2) account for material reference frame dependence under large 
deformations, and (3) satisfy first and second law thermodynamic constraints. The approach 
used here has been applied with success to model isotropic materials. 32 It is motivated by a 
focus on hypervelocity problems, where large deformation dynamics are of central interest, 
and by the kinematic form of the hybrid numerical method used in the present paper. As 
indicated in the later section on simulations, this material modeling approach provides an 
accurate description of hypervelocity impact effects in reinforced carbon-carbon. Potential 
application of the formulation to model impact in other composites, such as graphite-epoxy 
or Kevlar-epoxy, is of interest for future work. 

In the case of an orthotropic material, with distinct elastic moduli in tension and com- 
pression, the shear strain energy density per unit reference volume is 

= (1 - d) ftd £ fti [(1 + 7i) + (1 - 7i) sgn(Efr )} (ET? + 

2 = 1 
3 3 

(1 - d) ft, £ £ (1 - <%) ft, (EiP) 2 (18) 

4=1 j = 1 

where Sij is the Kronecker delta, /i 0 is a reference elastic modulus, and the parameters 
Hij = fiji are dimensionless constants. The parameters y* are the ratios of the elastic moduli 
in compression to those in tension, while the Ef™ are the components of the elastic shearing 
strain, expressed in a material reference frame. Note that this function is analytic, since 
a change in modulus from tension to compression occurs when the corresponding material 
strain component is zero. 

A plastic flow rule for an anisotropic, rate dependent material, which satisfies the afore- 
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mentioned isochoric plastic deformation constraint, may be obtained by extending a large 
strain Lagrangian formulation previously developed for use in hypervelocity impact applica- 
tions. 32,33 The flow rule is 

E,, = p^fSiM PSP ( 19 ) 

where A is a positive proportionality coefficient, S p is the effective stress, 


S" = m pT m t n t n pT s 


( 20 ) 


and S is the deviatoric stress tensor 

q _ 

<9E e 

The first two coefficients in the ffow rule impose the isochoric plastic deformation constraint, 
and are defined by 

= T= nT^ (c '" r + TC " ) (22) 

and 

gT = T— -tr(T)I (23) 

3 

for any symmetric second order tensor T. The third coefficient performs a component 
transformation from a material reference frame to a fixed global frame, and is defined by 

M T T = R t T R (24) 

for any symmetric second order tensor T. The last coefficient in the ffow rule defines an 
effective stress transformation, in a material reference frame, using 

M pT p = Q (25) 



for symmetric second order tensors P and Q, with component forms 

2 ota Pa 


Qii 


(1 + Pa) + (1 - Pu) sgn(Ef™ ) 


i = 1,2,3 


(26) 


146 


and 


Qij 


2 OLa Pj 


IJ * IJ 


(1 + fa) + (1 - fa) sgn(J - 1) 


3 


(27) 


The parameter a,j = atji is the ratio of a reference yield stress to the yield stress for the ijth 
stress component, while the parameter flij = (3ji is the ratio of the strength in compression 
to that in tension for the ijth stress component. 

The rate dependent, strain hardening, thermal softening yield function is 


/ 



(28) 


where Y is the yield stress 


Y = l (l-d)Y 0 (1 -k 9 h ) (1 + V £ p ) r 


i + c % ( — 

£os J 


(29) 


with Y 0 the reference yield stress, 77 a strain hardening coefficient, n a strain hardening 
exponent, ( a strain rate hardening coefficient, m a strain rate hardening exponent, e Q a 
reference strain rate, n a thermal softening coefficient, and 6 H the homologous temperature 

e-e n 


e H = ^^Y 

9 m ~ 9 0 


(30) 


where 9 0 and 9 m are reference and melt temperatures. 

In a numerical implementation, the aforementioned plastic flow rule is expressed in in- 
cremental form. That is the incremental plastic strain at each time step is computed using 
the incremental proportionality coefficient 

( ||S P || — Y \ 

AA = max 0 , 

V (1 — d) 2 Ho ) 

The shear damage variable (d) models the transition from an intact to a failed medium, 
evolving from an initial value of 0 to a final value of 1 over a fixed number of time steps 34 when 
any stipulated element failure criterion is satisfied. The simulations discussed in the next 
section incorporate accumulated plastic strain, melt temperature, and maximum compression 
failure criteria, although other criteria may be accommodated. 
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IMPACT SIMULATIONS 


The material model just described was applied in a series of three dimensional simulations 
of hypervelocity impacts on reinforced carbon-carbon. The simulations employed a hybrid 
particle-finite element method and the material properties listed in Tables 1 and 2. An 
initial set of simulations was used to validate the material model, compare results obtained 
using analytic (Mie Gruneisen) and tabular 35 (SESAME 3715) equations of state, and check 
numerical convergence of the simulation results. A second series of simulations was then 
performed to estimate orbital debris impact effects at velocities beyond the range of current 
experimental methods. 

The first set of eight simulations modeled NASA JSC experiments B1028 and B1040, 9 
which involved oblique impacts of aluminum spheres on reinforced carbon-carbon target 
plates at a velocity of seven kilometers per second. The target plates were 0.63 cm in 
thickness, including upper and lower surface coatings composed of silicon carbide, each 
0.08 cm in thickness. Table 3 lists the simulation parameters, including projectile diameter 
(D), impact velocity (u), impact obliquity (0, with zero degrees a normal impact), number 
of elements spanning the target thickness (N e ), and the equation of state used to model 
the aluminum projectile. Tabular equation of state data was not available for the target 
materials. 

Figures 1 through 4 show example plots for a simulation of experiment B1028. Figure 
1 shows the initial configuration while Figures 2 through 4 show the simulation results at 
50 microseconds after impact. The sectioned plot in Figure 4 depicts plate perforation 
and coating spall similar to that observed in the corresponding experiment. Table 3 lists 
simulation results for the diameter of the RCC perforation (D p ) and the average diameter 
of the target region over which the silicon carbide coating was removed (D c ). The results of 
the validation simulations suggest the following conclusions: 
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• the material model developed here can provide good estimates of both the RCC per- 
foration diameter and the extent of the spalled coating region, for oblique impacts at 
seven kilometers per second, 

• accurate estimates of the RCC perforation diameter require a mesh resolution sufficient 
to place 8 elements across the target plate, 

• accurate estimates of the diameter of the spalled coating region require a mesh resolu- 
tion sufficient to place 16 elements across the target plate, and 

• the simulation results are not sensitive to the choice of projectile equation of state. 

Table 4 shows the relative computational cost of simulations of experiments B1028 and 
B1040 run at three different mesh densities. As is well known, in three dimensional models 
the particle count increases with the cube of the increase in resolution, while the time step 
decreases linearly with the increase in resolution, so that the total computational cost of 
high resolution models is considerable. 

A second set of twelve simulations was performed to investigate orbital debris impact 
effects at velocities beyond the current experimental range. The simulations involved spher- 
ical aluminum projectiles, at three different projectile diameters, an impact obliquity of 30 
degrees and impact velocities of 7, 10, and 13 km/s. In the case of the largest projectile, 
simulations were performed using both an analytic and a tabular equation of state. The 
target assumed in these simulations was identical to that involved in the aforementioned ex- 
periments. In the target mesh 8 elements spanned the plate thickness, so that the resolution 
level was sufficient to estimate the diameter of the RCC perforations, but not the extent of 
the region of coating spall. Figure 5 shows simulation results for the diameters of the RCC 
perforations, as a function of projectile diameter, impact velocity, and projectile equation 
of state. In Figure 5, MG denotes the Mie Gruneisen analytic equation of state, while SES 
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denotes the SESAME tabular equation of state. The results of these simulations suggest the 
following conclusions, for the impact velocity range and impact obliquity considered: 

• perforation diameters increase with both projectile size and impact velocity, over the 
full range of the simulations, 

• for a fixed projectile size, perforation diameters increase with impact velocity at an 
approximately linear rate, 

• for a fixed impact velocity, perforation diameters increase with projectile size, but at 
a declining rate, and 

• the simulation results are not sensitive to the choice of projectile equation of state. 

Note that Figure 5 is not a ballistic limit plot; rather it plots perforation diameter versus 
impact velocity, so that the indicated trends are not unexpected. 

Although the preceding results are informative, they consider only a limited rage of 
projectile size and obliquity. Hence the scaling of the simulation results, as compared to the 
available experimental data, is of considerable interest. Figure 6 shows a plot of perforation 
diameter versus normal impact momentum for the 11 different projectile size and impact 
velocity combinations modeled in the present computational study, as well as corresponding 
data for 15 published experiments. 9,36 The experiments involved projectile diameters ranging 
from 0.039 to 0.628 cm, impact velocities ranging from 2.49 to 7.32 km/s, and impact 
obliquities ranging from 0 to 80 degrees. The simulations involved a more limited range 
of projectile diameters (0.123 to 0.360 cm) and obliquities (30 to 45 degrees), but a much 
higher range of impact velocities (7 to 13 km/s). All of the simulations and experiments 
of course involved the same target configuration. The data in Figure 6 suggests that the 
experimental and simulation results for the diameters of RCC perforations scale with normal 
impact momentum in a similar fashion. Although these results do not establish a universal 
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scaling relation for the problem of interest, they do suggest that the scaling of perforation size 
with normal impact momentum observed in experiments below 8 km/s may be extrapolated 
to much of the velocity range of interest in orbital debris impact applications. 

CONCLUSION 


The present paper has formulated an anisotropic, rate dependent material model for use in 
the simulation of hypervelocity impact problems. The material model was developed to study 
orbital debris impact effects on reinforced carbon-carbon materials, and has been validated in 
simulations of hypervelocity impact experiments conducted at 7 km/s. The validated model 
was applied to simulate impacts at velocities beyond the experimental range. The results 
indicate that momentum scaling analysis, used to correlate a wide range of experiments 
below 8 km/s, has application in predicting perforation diameters for reinforced carbon- 
carbon targets at velocities as high as 13 km/s. The ability of reinforced carbon-carbon 
to retain its strength at high temperatures suggests that accurate strength models of this 
material are important in simulations of impact effects over the entire orbital debris velocity 
range. 

Some conclusions relevant to future work are suggested: 

• additional high resolution simulations are needed in order to investigate the spallation 
of silicon carbide coating at velocities above the current experimental range, 

• additional mechanical properties testing is needed, at elevated temperatures and high 
strain rates, to support the development and validation of improved strength models 
for reinforced carbon-carbon, 

• additional equation of state research is needed, to provide tabular data applicable to 
reinforced carbon-carbon materials over a wide range of impact velocities, and 
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• the development of advanced thermal protection materials should in the future include 
experimental work aimed at detailed characterization of their mechanical as well as 
their thermal properties. 

The use of composite materials in spacecraft applications complicates both experimental 
and computational studies of impact effects. High cost materials with long fabrication lead 
times, such as reinforced carbon-carbon, and the limitations of current experimental impact 
techniques motivate the increased use of computer simulation in the design of spacecraft for 
micrometeoroid and orbital debris impact effects. 
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Table 1. Material properties 


Material property 

Aluminum 

Silicon Carbide 

RCC 

Reference bulk modulus ( Mbar ) 

0.784 

2.21 

0.0576 

Reference shear modulus (Mbar) 

0.271 

0.240 

0.0718 

Reference soundspeed ( cm /nsec -1 ) 

0.539 

0.829 

0.191 

Mie-Gruneisen gamma 

1.97 

0.95 

0.24 

Mie-Gruneisen slope 

1.34 

1.21 

1.33 

Reference density ( g cm~ 3 ) 

2.70 

3.21 

1.58 

Reference yield stress ( kbar ) 

2.90 

0.771 

0.771 

Specific heat ( bar cm 3 g~ x °K~ 1 ) 

8.84 

7.12 

7.12 

Strain hardening coefficient 

125 

10 

2 

Strain hardening exponent 

0.1 

1.0 

1.0 

Strain rate hardening coefficient 

0.0 

0.0 

0.1 

Strain rate hardening exponent 

0.0 

0.0 

1.0 

Reference strain rate (sec -1 ) 

0 

0 

0.01 

Thermal softening coefficient 

0.567 

0.0 

-1.0 

Melt temperature (°K) 

1,220 

3,840 

3,840 

Maximum compression 

100 

100 

100 

Plastic failure strain 

1.00 

0.10 

0.50 


Table 2. Material model parameters 


Parameters 

Silicon Carbide 

Reinforced carbon-carbon 

II 

bo 

II 

CO 

0.10 

1.00 

II 

to 

to 

II 

-£ 

CO 

CO 

II 

CO 

1.00 

1.00 

/R2 = /^23 

1.00 

0.50 

an = 022 = «33 = «13 

1.00 

1.00 

ai2 = Oi 23 

1.00 

3.73 

/Til = 7?22 = /?33 

2.00 

2.00 

/3l2 = /3l3 = /3 23 

2.00 

2.00 
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Table 3. Numerical results, simulations of NASA JSC experiments B1028 and B1040 


Test 

number 

D 

(cm) 

V 

(km/s) 

0 

(deg) 

N e 

Equation 
of state 

D p 

(cm) 

Error 

(%) 

D c 

(cm) 

Error 

(%) 

B1028 

0.628 

7.01 

45 

8 

Mie Gruneisen 

2.60 

10.3 

3.74 

15.0 





8 

SESAME 3715 

2.65 

8.6 

3.60 

18.2 





16 

Mie Gruneisen 

2.66 

8.3 

4.05 

8.0 





24 

Mie Gruneisen 

2.67 

7.9 

4.08 

7.3 

B1040 

0.478 

6.96 

30 

8 

Mie Gruneisen 

2.12 

3.6 

2.95 

21.3 





8 

SESAME 3715 

1.97 

10.5 

2.95 

21.3 





16 

Mie Gruneisen 

2.00 

10.0 

3.38 

9.9 





24 

Mie Gruneisen 

2.10 

4.5 

3.48 

7.2 


Table 4. Computer resource requirements, simulations of NASA JSC experiments 


Test 

number 

N e 

Total 

particles 

(million) 

Total 

elements 

(million) 

Number of 
processors 

Wall clock 
hours 

B1028 

8 

0.078 

0.036 

16 

14 


16 

0.572 

0.275 

64 

65 


24 

1.857 

0.905 

64 

340 

B1040 

8 

0.078 

0.036 

32 

5 


16 

0.572 

0.275 

64 

74 


24 

1.856 

0.905 

64 

347 
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Figure 1. Initial configuration, simulation of NASA JSC experiment B1028. 
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Figure 2. Particle-element plot of the simulation results at 50 microseconds after impact. 
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Figure 3. Element plot of the simulation results at 50 microseconds after impact. 
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Figure 4. Sectioned element plot of the simulation results at 50 microseconds after impact. 
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Figure 5. Simulation results for perforation diameter versus impact velocity. 
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Figure 6. Perforation diameter versus normal impact momentum, for hypervelocity impact 

in reinforced carbon-carbon. 
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ABSTRACT 


A series of three dimensional simulations has been performed to investigate analytically the effect of 
insulating foam impacts on ceramic tile and reinforced carbon-carbon components of the Space Shuttle 
thermal protection system. The simulations employed a hybrid particle-finite element method and a parallel 
code developed for use in spacecraft design applications. The conclusions suggested by the numerical study 
are in general consistent with experiment. The results emphasize the need for additional material testing 
work on the dynamic mechanical response of thermal protection system materials, and additional impact 
experiments for use in validating computational models of impact effects. 

NOMENCLATURE 

c particle center of mass position vector 
d element shear damage variable 

D element normal damage variable 

E p element plastic strain tensor 

f particle damping force 

m particle mass 

p particle translational momentum vector 

S particle entropy 

V thermomechanical potential energy function 

INTRODUCTION 

the Columbia Accident Investigation Board 1 concluded that the effects of foam impact on 
edge of the Space Shuttle was the most likely cause for the loss of the Orbiter Columbia. 
Strong evidence in support of this conclusion is provided by a recent series of impact experiments conducted 
at Southwest Research Institute (SwRI) by a NASA-SwRI-industry team . 2 The current consensus regarding 
the cause of the accident was not present in the early stages of the investigation, since little experimental data 
relevant to the accident conditions was available, and since significant lead times were required to prepare 
and conduct the necessary impact experiments. 


The report of 
the wing leading 
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Soon after the loss of Columbia an impact analysis team was assembled 3 whose purpose was to investigate 
analytically the effects of foam impact on components of the Space Shuttle thermal protection system, and 
to support the conduct of experiments designed to duplicate the impact events observed during launch 
of the vehicle. This group included NASA, industry, national laboratory, and university participation and 
employed a variety of numerical methods 4 and computer codes 5 to simulate the impact events of interest. The 
present paper describes simulations performed using a hybrid particle-finite element method 6 and a parallel 
computer code 7 to model the impact of foam blocks on both ceramic tile and reinforced carbon-carbon 
(RCC) components of the Space Shuttle thermal protection system. The simulations described here were 
performed in advance of the aforementioned experiments and employed the best available material property 
data for foam, tile, and RCC. The conclusions suggested by the simulations are in general consistent with the 
results of later experiments, although additional material testing, material modeling, and simulation work 
is needed to develop a validated computational approach to impact damage assessment for future Space 
Shuttle applications. 

The sections which follow describe the numerical method used in the simulations, the structural and 
material models assumed for the foam projectiles and the ceramic tile or RCC targets, the computational 
costs of the simulations, and the results of the numerical study, including suggestions for future research. 

NUMERICAL METHOD 

In recent research focused on the design of orbital debris shielding, a new numerical method and parallel 
computer code have been developed for use in spacecraft design applications. This hybrid numerical method 
employs in tandem nondeforming Lagrangian particles and large strain finite element kinematics , 8 to simulate 
impact problems involving shock loading, large deformation plasticity, and complex fragmentation dynamics. 
The method has been implemented in a three dimensional code and validated by comparison with published 
experiments at impact velocities ranging from one to eleven kilometers per second . 9 

The hybrid method combines the general contact-impact modeling capabilities of particle methods with 
a true Lagrangian description of strength effects, the latter offered by finite element techniques. It avoids 
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the tensile instability problems which have hindered the effective use of some particle techniques, as well as 


the mass and energy discard normally associated with Lagrangian finite element models of material failure. 
No particle to element mapping is required, since both particles and elements are used throughout the 
calculation, but to represent distinct physical effects. The particles model all inertia and all contact-impact 
as well as volumetric thermomechanical response in compressed states, while the elements model tension and 
elastic-plastic shear. Material failure is represented by the loss of element cohesion, after which particles not 
associated with any intact element are free to flow under general contact-impact loads. 

In the case of spherical particles, the state space model for the particle-element system 8 consists of 
evolution equations for the particle translational momenta and center of mass position vectors (c^) 


,(0 = _ f (*) 


P w = - 


<9cW 


c W = toW -1 p (i) 


(1) 


augmented by evolution equations for the internal state variables 



= S {i) (p« , c« , S (i) , d {j) , D & , E p «) 

(2) 

d u) 

= d {j) (p W , c W , S (i) , d U) , D & , E p(i) ) 

(3) 

i)U) 

= Z)k') ( p (d ; c (i) , S {i) , d U) , D (i) , E p(j) ) 

(4) 


= E p W (p (i) , c (i) , S (l) , d U) , D ^ , E pw ) 

(5) 


where f W is a damping force, m W is a particle mass, is a particle entropy, d^ and are element 
shear and normal damage variables, E p ^'^ is a plastic strain tensor, i is a particle index, j is an element 
index, and V is a thermomechanical potential 


V = V(c (i) ,S (i) ,d U) ,D u) ,E p(j) ) 


( 6 ) 


The specific functional forms of the thermomechanical potential and the internal state evolution equations 
depend upon the constitutive assumptions as well as the adopted interpolations for the density and displace- 
ment fields. The present work investigated for the first time the application of this method to a relatively low 
velocity impact regime, in problems which nonetheless involved complex contact-impact, material failure, 
and fragmentation phenomena difficult to simulate using structural finite element codes. 
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MATERIAL MODELS 


The simulations described here employed simple material models for the foam, tile, felt strain isolation 
pad (SIP), and RCC, with material properties estimated using the available experimental data. All materials 
were assumed to be isotropic elastic-perfectly plastic, 10 with an accumulated a plastic strain criterion applied 
to initiate element failure. The available material data base may be summarized as follows. In support of 
the Columbia accident investigation, Glenn Research Center 11 and Sandia National Laboratories (SNL) 12 
performed mechanical property tests on foam, tile, and reinforced carbon-carbon. Mechanical property 
tests previously performed on SIP and on SIP-tile combinations are described by Sawyer 13 and Cooper 
and Sawyer 14 respectively. The relevant thermal properties for polyurethane, tile, SIP, and carbon-carbon 
materials are provided by Oertel, 15 Banas et al., 16 Myers et al., 17 Ohlhorst et al., 18 and the commercial 
literature. 19 

Table 1 lists estimated properties for the materials of interest. These values were used (except in the case 
of the SIP) to perform the simulations described in this paper. The present analysis adopted a yield strength 
for the tile corresponding to the lowest experimental measurements of Lu et al. 12 In the case of the RCC, 
the present analysis assumed a yield strength equal to the bending strength measured by Lu et al., 12 since 
the failure mode for the RCC panels was expected to be flexure of the panel surface under the foam impact 
load. At the time the present analysis was performed, the data available to the authors on SIP properties 
was very limited. As a result the SIP density was underestimated by a factor of 2.3, and the single layer of 
SIP elements used in the numerical model was assigned the same stiffness properties as the tile. The effect 
of the underestimating the SIP density was to slightly underestimate the target areal density, and in the 
present analysis is not considered to be significant. The experiments of Cooper and Sawyer 14 suggest that 
the stiff SIP elements used here would tend to overestimate the tile damage produced by the foam impact 
load. 

More general models of the dynamic mechanical deformation and failure of the foam, tile, felt, and 
reinforced carbon-carbon materials are needed. The authors and others are currently engaged in work to 
develop improved material models, with anisotropy and strain rate dependence significant in some if not all 
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of these materials. 


TILE IMPACT MODEL 

Analysis of launch videos, supplemented by computational fluid dynamics studies, suggested that ceramic 
tiles located on the lower surface of Columbia might have been struck by a block of insulating foam shed from 
the vehicle’s external tank. A series of experiments was therefore planned to measure the impact damage 
produced by highly oblique foam block impacts on tile arrays similar to those covering the lower surface of 
the Orbiter wing. Simulations were run in advance of these experiments, to estimate the impact damage. 
The simulation parameters are listed in Table 2. 

In each simulation the target model was composed of a 2 x 4 foot array of tiles, each tile having an areal 
extent of 6 x 6 inches. The uniform tile thickness matched that in the suspected impact areas, over the main 
landing gear door and in the nearby wing acreage. In the simulations the tile array was supported by an 
aluminum plate, whose lower surface was fixed along a circumferential edge strip, the latter with a width 
of one inch. The strain isolation pad (SIP) which separates the tile and the aluminum wing structure was 
modeled with a single layer of finite elements, however similar (gap filler) material often interposed between 
the individual tiles was not modeled. The foam projectile was modeled as a homogeneous hexahedral block. 
The dimensions, obliquity, and orientation of the foam block at impact were varied between simulations, 
due to uncertainties in the interpretation of the launch videos, a dependence of the impact obliquity on the 
vehicle impact location, and a desire to investigate the effect of projectile orientation (roll angle) on impact 
damage. 

Computer resource requirements and some limitations of the research code and preprocessor used here 
made it necessary to introduce certain geometric approximations. Since available commercial preprocessors 
do not generate hybrid particle-finite element models, a special preprocessor was employed. The latter code 
generates solid models composed of uniform hexahedra, and associated ellipsoidal particles, so that an element 
and particle deletion process was used to introduce the gaps between the tiles. As a result the width of these 
gaps was overestimated. Combined with the aforementioned neglect of gap fillers, the assumed geometric 
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model approximates conservatively the structural strength of the actual tile array. A second approximation 
was introduced in modeling the individual tiles, whose external surfaces are coated in a borosilicate layer, to a 
depth approximately five percent of the tile thickness. Computer resource requirements precluded modeling 
of features with such small dimensions, so the tiles were taken to be monolithic with material properties 
derived from the published strength and stiffness properties of an individual tile. 

TILE IMPACT SIMULATIONS 

The four tile impact simulations were performed on systems operated by NASA Ames Research Center, 
and required between three and five wall clock days on 128 or 192 processors of an SGI Origin. The models 
were composed of over one million particles, with the simulations extending over five or six milliseconds of 
physical time. 

The first two simulations differed only with respect to projectile orientation (roll angle), and modeled the 
impact of a 1.06 pound block of foam, at a velocity of 700 feet per second, on a tile array similar to those 
which cover the main landing gear doors. Impact obliquity was five degrees. The simulations showed 12.0 
cubic inches of material (0.178 tile volumes) eroded by the long edge impact and 19.6 cubic inches of material 
(0.290 tile volumes) eroded by the short edge impact. Figures 1 and 2 show the simulation results, including 
views of the predicted tile erosion. In both these simulations the maximum predicted depth of penetration 
was approximately one half of the tile thickness. Since the short edge impact appeared to be more damaging 
under the postulated impact conditions, subsequent simulations (and later experiments) involved foam blocks 
rotated so as to strike the tile array along the projectile’s short edge. Experiments which approximately 
correspond to these two simulations were later performed by Kerr et al. 2 The first experiment produced 
three craters with a total volume of approximately 0.1 cubic inches, although much of this damage may have 
been caused by the unintended impact of a Mylar burst disc used in the compressed air gun which launched 
the projectile. The second experiment produced no impact craters. Since the aforementioned numerical 
modeling assumptions minimized both tile strength and tile lateral support while maximizing the stiffness 
of the SIP layer, the tile impact damage was overestimated. However the eroded volume error was less than 
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one third of one tile volume, with the damage distributed among several tiles. 

The third and fourth simulations considered somewhat more severe impact conditions. The third case 
assumed a 2.24 pound projectile and a slightly reduced tile thickness (tile thickness varies with position 
over the lower surface of the Orbiter). The result was an increase in eroded material, to 48.2 cubic inches 
(0.789 tile volumes) and an increase in the maximum depth of penetration, to three quarters of the tile 
thickness. These simulation results are depicted in Figure 3. Finally the fourth case considered a lower 
mass projectile (1.53 pounds), but a slightly higher impact velocity (720 fps) and a more direct impact, with 
the impact obliquity taken to be thirteen degrees. The tile thickness in the target was increased in order 
to represent a wing acreage area away from the main landing gear door. The simulation results, shown in 
Figure 4, predicted erosion of 70.6 cubic inches of material (0.785 tile volumes) and a maximum depth of 
penetration to the level of the SIP, in one small area. Experiments which approximately correspond to these 
two simulations were later performed by Kerr et al. 2 The first experiment produced no impact craters. The 
second experiment produced four craters, each with an areal extent of less than 1.0 square inches (depth not 
provided) . It appears that the aforementioned numerical modeling assumptions again caused the tile impact 
damage to be overestimated. The eroded volume error was less than one tile volume, with the damage 
distributed among multiple tiles. 

In summary the pretest simulations predicted in the worst case the removal of less than one tile volume of 
ceramic material, under modeling assumptions which conservatively approximated tile strength properties, 
the lateral support provided in the tile gap region, and the compliance of the SIP layer. Subsequent testing 
showed that none of the impact configurations considered here produced significant damage to the target 
tile array. 


RCC IMPACT MODEL 

As in the case of the underwing tiles, analysis of launch videos and complimentary computational fluid 
dynamics work suggested the possibility that Columbia’s wing leading edge was subjected to a highly oblique 
foam block impact. A series of experiments on reinforced carbon-carbon panels was therefore planned to 
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investigate the effects of such impacts, for panel geometries representative of the leading edge region most 
likely involved. Prior to these experiments, two simulations were run to estimate the impact damage. The 
simulation parameters are listed in Table 3. 

The target model used in the simulations represented the geometry of wing leading edge panel number 
six. The limitations of the preprocessor used here again led to certain approximations. A profile for the 
model cross section was obtained by fitting coordinate data extracted from a CAD model of the actual 
panel, and assuming a constant RCC wall thickness. This cross section was then extended an axial distance 
equal to the total panel length, with stiffening ribs added at both ends, similar to those found on the actual 
part. The upper and lower edges of the panel were held fixed in the simulations. This target model was 
considered to be generally representative of the strength and stiffness of the actual structure. The simulations 
assumed that in the RCC elements failure would occur at a plastic strain of 0.01, a relatively brittle failure 
criterion. The thin silicon carbide coating present on the actual part was not modeled, again due to the 
high computational cost of simulations which resolve very small scale features. As discussed in a preceding 
section, the particle-element preprocessor used here produced models composed of uniform hexahedra, so 
that the curved surface of the RCC panel model was represented with a stairstep geometry. 

The starting conditions for the two simulations differed only with respect to projectile orientation (roll 
angle), one objective of the analysis being to determine the relative severity of impact damage caused by edge 
and corner impacts. The specified impact point was located a distance of 18.1 inches from the panel edge, 
measured along the panel arc, and the impact obliquity (14.6 degrees) was specified as the angle between 
the target surface normal at the impact point and the projectile velocity vector, the latter aligned with the 
long axis of the foam block. The pitch, roll, and yaw of the projectile were computed so as to match these 
specifications. 


RCC IMPACT SIMULATIONS 

The RCC panel impact simulations were performed on SGI Origin systems at NASA Ames Research 
Center, and required between three and four wall clock days on 128 or 256 processors. The models were 
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composed of approximately two million particles, with the simulations extending over no more than two 
milliseconds of physical time. The simulation results are depicted in Figure 5 (corner impact case) and 
Figure 6 (edge impact case). The first simulation modeled a corner impact and resulted in failure of the 
panel, with a crack approximately six inches in length developing along the panel surface, normal to the 
leading edge stagnation line. The second simulation modeled an edge impact and showed greater panel 
damage, in this case multiple cracks, the largest extending along half the length of the panel and aligned 
parallel to the stagnation line. 

Although the scope of the RCC impact simulation work described here was limited, the results indicate 
failure of the panel under the postulated foam impact loads. These results are in general consistent with 
later experiments conducted on Space Shuttle wing leading edge panels. Experiments which approximately 
correspond to these two simulations were performed by Kerr et al. 2 In the first experiment, a corner impact 
test conducted on a panel six target at a twenty-one degree obliquity resulted in a crack of length 5.5 inches 
located at the panel edge and oriented parallel to the stagnation line. In the second experiment, an edge 
impact test conducted on a panel eight target at a twenty-five degree obliquity resulted in gross failure of the 
panel surface, producing a 17 inch by 16 inch hole in the panel surface. Comparison of the experiments and 
simulations is complicated by differences in target geometry and impact obliquity. Since the experiments 
involved higher impact obliquities, they would be expected to produce more damage than is depicted in the 
simulations. The first simulation showed a crack similar in size to that observed in the first experiment, 
although the predicted location and orientation were not correct. The second simulation showed large cracks 
in the panel surface, at an impact obliquity ten degrees less than that which resulted in gross panel failure in 
the second experiment. In general the numerical modeling work, which incorporated best case assumptions 
with regard to RCC strength and ductility, appears to provide good estimates of panel impact damage. 

CONCLUSION 

The present paper has described a series of pre-test simulations performed to estimate damage produced 
by external foam strikes on thermal protection system components of the Space Shuttle. The simulations 
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employed a hybrid particle-finite element technique and a parallel computer code developed for use in space- 
craft design applications. The simulation results are in general consistent with experimental results available 
for this class of problems, and indicate that the numerical method used here is suitable for application in 
a relatively low velocity regime. The application of this numerical technique to future impact problems 
would be facilitated by further methods and interface development work, aimed at accommodating complex 
structural geometries described by a standard CAD data base or a commercial finite element preprocessor. 

Several conclusions specific to the operation of the Space Shuttle and the design of future aerospace planes 
are suggested: (1) additional material testing and constitutive modeling research describing the deformation 
and failure of thermal protection system materials is needed, (2) numerical methods and code development 
work is needed to provide a validated computer simulation capability for impact damage assessment, (3) 
additional, higher resolution simulations should be performed, to investigate the effects of any simplifying 
assumptions made in the areas of material modeling and structural geometry, and (4) additional impact 
testing should be conducted, over a wider range of impact conditions, to validate proposed computational 
analysis techniques. 

Research in the suggested areas is already in progress. The authors are currently engaged in work which 
will allow the hybrid particle-finite element technique used here to model impact on any structural geometry 
described by a general hexahedral finite element mesh. 
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Table 1. Material properties 


Material property 

Foam 

Tile 

SIP 

RCC 

Young’s modulus (psi) 

1,360 

9,022 

220 

2.21x 10® 

Shear modulus (psi) 

529 

3,510 

110 

1.04x10® 

Reference density (lb /ft 3 ) 

2.15 

9.00 

12.3 

98.6 

Yield stress (psi) 

42.2 

23.9 

27.4 

14.0x10 s 

Specific heat (Btu/lbm-degree F) 

0.454 

0.151 

0.315 

0.171 

Thermal expansion coefficient (1 /degree F) 

0 

2.25x 10~ 7 

l.OOxlO -5 

7.28x 10~ 7 

Plastic failure strain 

1.0 

1.0 

1.0 

0.01 


Table 2. Parameters of the tile impact simulations 


Parameter 

MLGD 1 

MLGD 2 

MLGD 3 

Wing Acreage 

Projectile velocity (fps) 

700 

700 

700 

720 

Impact obliquity (degrees) 

5 

5 

5 

13 

Projectile roll (degrees) 

90 

0 

0 

0 

Projectile cross section (in) 

3.5x11.5 

3.5x11.5 

5.5x11.5 

5.5x11.5 

Projectile length (in) 

21.25 

21.25 

28.5 

19.0 

Tile thickness (in) 

1.875 

1.875 

1.700 

2.450 

Aluminum plate thickness (in) 

0.1875 

0.1875 

0.2689 

0.2720 

Simulation time (milliseconds) 

5.0 

5.0 

6.0 

5.0 

Number of particles (millions) 

1.10 

1.10 

1.42 

1.48 

Number of processors (SGI Origin) 

128 

128 

192 

192 

Wall clock time (hours) 

76.5 

90.1 

127 

102 


Table 3. Parameters of the wing leading edge impact simulations 


Parameter 

Comer impact case 

Edge impact case 

Projectile velocity (fps) 
Projectile dimensions (in) 
Projectile roll, pitch, yaw (degrees) 
Impact obliquity (degrees) 
Panel dimensions (in) 

Panel thickness (in) 
Simulation time (milliseconds) 
Number of particles (millions) 
Number of processors (SGI Origin) 
Wall clock time (hours) 

775 

5.5x11.5x22.8 
0.0, 17.5, 6.32 
14.6 

20.5x38.4x21.3 

0.25 

1.635 

1.90 

128 

96 

775 

5.5x11.5x22.8 
30.0, 17.5, 6.32 
14.6 

20.5x38.4x21.3 

0.25 

2.000 

1.90 

256 

74 
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Figure 1. Simulation results for case MLGD 1. 
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Figure 2. Simulation results for case MLGD 2. 
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Figure 3. Simulation results for case MLGD 3. 



181 






Figure 4. Simulation results for the wing acreage impact case. 
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Figure 5. Simulation results for the RCC corner impact case. 
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Figure 6. Simulation results for the RCC edge impact case. 
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CONCLUSIONS 


The design of manned spacecraft for future space exploration missions will require con- 
sideration of micrometeoroid and orbital debris impact effects. The debris environment in 
low earth orbit presents a significant hazard, and has mandated the development of hy- 
pervelocity impact shielding for the International Space Station. Although exposure times 
for the Space Shuttle are much less than that for the space station, the shuttle routinely 
suffers limited orbital debris impact damage. Hence next generation spacecraft intended to 
operate, even for limited periods, in low earth orbit must be designed with the orbital de- 
bris impact threat in mind. For operations beyond earth orbit, the principal impact threat 
is due to micrometeoroids. Micrometeoroid impacts typically involve lower particle masses 
and densities, but higher impact velocities, than those associated with man-made debris in 
low earth orbit. As a result, shielding designed for low earth orbit projectiles may not be 
optimal to address the significant micrometeoroid hazard. Although considerable previous 
work has focused on the hypervelocity impact shielding problem, the development of next 
generation spacecraft is likely to require significant new experimental and computational re- 
search efforts. There are two principal reasons. The first is that the current knowledge base 
is focused heavily on projectiles, shielding, and spacecraft structures composed of metals. 
The second is that current experimental methods are not able to address the full impact 
velocity and kinetic energy range of interest. Although previous experimental and compu- 
tational studies of debris shielding have involved composites, work to date has served to 
highlight the increased complexity and cost of experimental and computational impact work 
involving advanced materials. The research described here on hypervelocity impact effects 
in reinforced carbon-carbon illustrates a design methodology likely to apply in future devel- 
opment of manned vehicles for space exploration missions. Coordinated experimental and 
computational efforts will likely be required to address orbital debris and micrometeoroid 
related design requirements for new space exploration systems. 
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APPENDIX 1 


Simulation Data for Reinforced Carbon-Carbon 


Simulation 

number 

D 

(cm) 

V 

(krn/s) 

0 

(deg) 

N e 

Equation 
of state 

D p 

(cm) 

Error 

(%) 

D c 

(cm) 

Error 

(%) 

200 

0.628 

7.01 

45 

8 

Mie Gruneisen 

2.60 

10.3 

3.74 

15.0 

219 




8 

SESAME 3715 

2.65 

8.6 

3.60 

18.2 

211 




16 

Mie Gruneisen 

2.66 

8.3 

4.05 

8.0 

218 




24 

Mie Gruneisen 

2.67 

7.9 

4.08 

7.3 

201 

0.478 

6.96 

30 

8 

Mie Gruneisen 

2.12 

3.6 

2.95 

21.3 

220 




8 

SESAME 3715 

1.97 

10.5 

2.95 

21.3 

212 




16 

Mie Gruneisen 

2.00 

10.0 

3.38 

9.9 

221 




24 

Mie Gruneisen 

2.10 

4.5 

3.48 

7.2 

191 

0.123 

7 

30 

8 

Mie Gruneisen 

0.159 

na 

nr 

na 

192 


10 


8 

Mie Gruneisen 

0.390 

na 

nr 

na 

193 


13 


8 

Mie Gruneisen 

0.536 

na 

nr 

na 

188 

0.240 

7 

30 

8 

Mie Gruneisen 

0.95 

na 

nr 

na 

189 


10 


8 

Mie Gruneisen 

1.21 

na 

nr 

na 

190 


13 


8 

Mie Gruneisen 

1.54 

na 

nr 

na 

185 

0.360 

7 

30 

8 

Mie Gruneisen 

1.41 

na 

nr 

na 

186 


10 


8 

Mie Gruneisen 

1.84 

na 

nr 

na 

187 


13 


8 

Mie Gruneisen 

2.22 

na 

nr 

na 

224 

0.360 

7 

30 

8 

SESAME 3715 

1.50 

na 

2.00 

na 

223 


10 


8 

SESAME 3715 

1.83 

na 

2.75 

na 

222 


13 


8 

SESAME 3715 

2.33 

na 

3.17 

na 


na = not available (no corresponding experiment) 
nr = not recorded (accurate coating spall results require N e = 16) 

D = projectile diameter (aluminum sphere) 
v = impact velocity; 0 = impact obliquity (normal impact is zero degrees) 
N e = number of elements across the target plate 
D p = perforation diameter; D c = diameter of the spalled coating region 
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